Open Case Studies : Predicting Annual Air Pollution

Disclaimer: The purpose of the Open Case Studies project is to demonstrate the use of various data science methods, tools, and software in the context of messy, real-world data. A given case study does not cover all aspects of the research process, is not claiming to be the most appropriate way to analyze a given data set, and should not be used in the context of making policy decisions without external consultation from scientific experts.

Motivation

A variety of different sources contribute different types of pollutants to what we call air pollution. Some sources are natural while others are anthropogenic (human derived):

[source]

Major types of air pollutants

1. Gaseous - Carbon Monoxide (CO), Ozone (O₃), Nitrogen Oxides(NO, NO₂), Sulpher Dioxide (SO₂)
2. Particulate - small liquids and solids suspended in the air (includes lead- can include certain types of dust)
3. Dust - small solids (larger than particulates) that can be suspended in the air for some time but eventually settle
4. Biological - pollen, bacteria, viruses, mold spores

See [here])http://www.redlogenv.com/worker-safety/part-1-dust-and-particulate-matter) for more detail on the types of pollutants in the air.

Particulate pollution

Air pollution particulates are generally described by their size.

There are 3 major categories:

1. Large Coarse Particulate Mater - has diameter of >10 micrometers (10 µm)

2. Coarse Particulate Mater (called PM_10-2.5) - has diameter of between 2.5 µm and 10 µm

3. Fine Particulate Mater (called PM_2.5) - has diameter of < 2.5 µm

PM₁₀ includes any particulate mater <10 µm (both coarse and fine particulate mater)

Here you can see how these sizes compare with a human hair:

[source]

The following plot and table show the relative sizes of these different pollutants in micrometers(µm):

[source]

[source]

This table shows how deeply some of the smaller fine particles can penetrate within the human body:

[source]

Negative Impact of Particulate Exposure on Health

Exposure to air pollution is associated with higher rates of mortality in older adults and is known to be a risk factor for many diseases and conditions including but not limited to:

1. Asthma - fine particle exposure (PM_2.5) was found to be associated with higher rates of asthma in children
2. Inflammation in type 1 diabetes - fine particle exposure (PM_2.5) from traffic-related air pollution was associated with increased measures of inflammatory markers in youths with type 1 diabetes
3. Lung function and emphysema - higher concentrations of ozone (O₃), nitrogen oxides (NO_x), black carbon, and fine particle exposure PM_2.5 , at study baseline were significantly associated with greater increases in percent emphysema per 10 years
4. Low birthweight - fine particle exposure(PM_2.5) was associated with lower birth weight in full-term live births
5. Viral Infection - higher rates of infection and increased severity of infection are associated with higher exposures to pollution levels including fine particle exposure (PM_2.5)

See this review article for more information about sources of air pollution and the influence of air pollution on health.

Sparse Monitoring is Problematic for Public Health

Historically epidemiological studies would assess the influence of air pollution on health outcomes by relying on a number of monitors located around the country. However as can be seen in the following figure, these monitors remain to be relatively sparse in certain regions of the country. Furthermore, dramatic differences in pollution rates can be seen even within the same city.

[source]

This lack of granularity in air pollution monitoring has hindered our ability to discern the full impact of air pollution on health and to identify at-risk locations.

Machine Learning Offers a Solution

An article published in the Environmental Health journal dealt with this issue by using data about population density, road density, among other features to model or predict air pollution levels at a more localized scale using machine learning methods.

Yanosky, J. D. et al. Spatio-temporal modeling of particulate air pollution in the conterminous United States using geographic and meteorological predictors. Environ Health 13, 63 (2014).

The authors of this article state that:

“Exposure to atmospheric particulate matter (PM) remains an important public health concern, although it remains difficult to quantify accurately across large geographic areas with sufficiently high spatial resolution. Recent epidemiologic analyses have demonstrated the importance of spatially- and temporally-resolved exposure estimates, which show larger PM-mediated health effects as compared to nearest monitor or county-specific ambient concentrations.”

The article above explains that machine learning methods can be used to predict air pollution levels when traditional monitoring systems are not available in a particular area or when there is not enough spatial granularity with current monitoring systems. We will use similar methods to predict annual air pollution levels spatially within the US.

Main Questions

Our main question:

1. Can we predict annual average air pollution concentrations at the granularity of zip code regional levels using predictors such as data about population density, urbanization, road density, as well as, satellite pollution data and chemical modeling data?

Learning Objectives

In this case study, we will walk you through importing data from CSV files and performing machine learning methods to predict our outcome variable of interest (in this case annual fine particle air pollution estimates). We will especially focus on using packages and functions from the Tidyverse, and more specifically the tidymodels package/ecosystem primarily developed and maintained by Max Kuhn and Davis Vaughan. This package loads more modeling related packages like rsample, recipes, parsnip, yardstick, and dials. We will also briefly cover the workflows and tune packages. The tidyverse is a library of packages created by RStudio. While some students may be familiar with previous R programming packages, these packages make data science in R especially efficient.

We will begin by loading the packages that we will need:

library(here)
library(readr)
library(dplyr)
library(skimr)
library(summarytools)
library(magrittr)
library(corrplot)
library(RColorBrewer)
library(ggcorrplot)
library(GGally)
library(tidymodels)# broom, dials, infer, parsnip, purrr, recipes, rsample, tibble, yardstick
library(workflows)
library(vip)
library(tune)
library(randomForest)
library(ggplot2)
library(stringr)
library(tidyr)



# library(stringr)
# library(purrr)
# library(tibble)
# library(tidyr)
# library(ggplot2)
# library(ggpubr)
# library(forcats)
# library(lmerTest)
# library(car)
# library(ggiraph)
# library(ggforce)
# library(viridis)
# library(cowplot)

Package	Use
here	to easily load and save data
readr	to import the CSV file data
dplyr	to view/arrange/filter/select/compare specific subsets of the data
skimr	to get an overview of data
summarytools	to get an overview of data in a different style
magrittr	to use the %<>% pipping operator
corrplot	to make large correlation plots
ggcorrplot	also to make large correlation plots
GGally	to make smaller correlation plots
rsample	to split the data into testing and training sets and to split the training set for cross-validation
recipes	to pre-process data for modeling in a tidy and reproducible way and to extract pre-processed data (major functions are recipe() , prep() and various transformation step_*() functions, as well as juice() - extracts final preprocessed training data and bake() - applies recipe steps to testing data). See here for more info.
parsnip	an interface to create models (major functions are fit(), set_engine())
yardstick	to evaluate the performance of models

stringr | to manipulate the text within the PDF of the data

purrr | to perform functions on all columns of a tibble tibble | to create data objects that we can manipulate with dplyr/stringr/tidyr/purrr tidyr | to separate data within a column into multiple columns ggplot2 | to make visualizations with multiple layers ggpubr | to easily add regression line equations to plots forcats | to change details about factors (categorical variables) lmerTest| to perform linear mixed model testing car| to perform Levene’s Test of Homogeneity of Variances ggiraph| to make plots interactive ggforce| to modify facets in plots viridis| to plot in color palette that is easily interpreted by colorblind individuals cowplot | to allow plots to be combined ___

The first time we use a function, we will use the :: to indicate which package we are using. Unless we have overlapping function names, this is not necessary, but we will include it here to be informative about where the functions we will use come from.

Context

The State of Global Air is a report released every year to communicate the impact of air pollution on public health.

The State of Global Air 2019 report which uses data from 2017 stated that:

Air pollution is the fifth leading risk factor for mortality worldwide. It is responsible for more deaths than many better-known risk factors such as malnutrition, alcohol use, and physical inactivity. Each year, more people die from air pollution–related disease than from road traffic injuries or malaria.

The report also stated that:

In 2017, air pollution is estimated to have contributed to close to 5 million deaths globally — nearly 1 in every 10 deaths.

##### [source]

The State of Global Air 2018 report using data from 2016 which separated different types of air pollution, found that particulate pollution was particularly associated with mortality.

The 2019 report shows that the highest levels of fine particulate pollution occurs in Africa and Asia and that:

More than 90% of people worldwide live in areas exceeding the World Health Organization (WHO) Guideline for healthy air. More than half live in areas that do not even meet WHO’s least-stringent air quality target.

Looking at the US specifically, air pollution levels are generally improving. The US Environmental Protection Agency (EPA) also releases a report about air pollution levels called Our Nation’s Air.

[source]

However, air pollution continues to contribute to health risk for Americans, in particular in regions with higher than national average rates of pollution that actually at time exceed the world health organization’s recommended level. Thus it is necessary to obtain high spatial granularity in estimates of air pollution in order to identify locations where populations are experiencing harmful levels of exposure.

You can see that current air quality conditions at this website and you will notice variation across different cities.

Here were the conditions in Topeka Kansas when this was written:

It reports particulate values using what is called the Air Quality Index scale (AQI), this calculator indicates that 114 AQI is equivalent to 40.7 ug/m³ and is considered unhealthy for sensitive individuals. Thus some areas very much exceed the World Health Organization (WHO) annual exposure guideline (10 ug/m³) at certain times and this may adversely affect the health of people living in these locations.

Furthermore, adverse health effects have been associated with populations experiencing higher pollution exposure despite the levels being below suggested guidelines. Secondly, it appears that the composition of the particulate mater and the influence of other demographic factors may make specific populations more at risk for adverse health effects due to air pollution. See this article for more details.

The monitor data that we will use in this case study comes from a system of monitors in which roughly 90% are located within cities. Thus there is an equity issue in terms of capturing the air pollution levels of more rural areas. Therefore, to get a better sense of the pollution exposures for the individuals living in these areas, methods like machine learning can be very useful to estimate air pollution levels in areas with little to no monitoring.

Indeed, machine learning methods are in fact used to be able to estimate air pollution in these low monitoring areas so that we can make a map like this where we have annual estimates for all of the contiguous US:

[source]

This is what we aim to achieve in this case study.

Limitations

There are some important considerations regarding this data analysis to keep in mind:

1. The data in our analysis does not include information about the composition of particulate mater. Different types of particulates may be more benign or deleterious for health outcomes.

2. Outdoor pollution levels are not necessarily an indication of of individual exposures. People spend differing amounts of time indoors and outdoors and are exposed to different pollution levels indoors. People are now developing personal monitoring systems to track air pollution levels on the personal level.

Our analysis will use annual mean estimates, however pollution levels can vary greatly by season, day and even hour. There are data sources that have finer levels of temporal data, however we are interested in long term exposures, as these appear to be the most influential for health outcomes, so we chose to use annual level data.

What are the data?

In Machine Learning for prediction, there are two main types of variables:

1. Outcome variable
2. Predictor variables

The outcome variable is what are trying to predict. In building our model we actually have the outcome variable data, but we want to see how well our predictor variables can explain the variation in our outcome data. This gives us a sense of how well we can use the predictor variable data to predict our outcome variable levels when we in fact do not have data about the outcome.

As a simpler example, imagine that we have data about the sales and characteristics of cars from last year and we want to predict which cars might sell well this year. We do not have the sales data yet for this year, but we do know the characteristics of our cars for this year. We can use a model of the characteristics that explained sales last year to estimate what cars might sell well this year. In this case, our outcome variable is the sale performance of the cars, while the different characteristics of the cars make up our predictor or explanatory variables.

In this case study, we will evaluate air pollution monitor data of fine particulate mater (PM_2.5) in the contiguous US from 2008, as well as data about population density, road density, urbanization levels, and NASA satellite data to develop models to predict localized air pollution levels.

The monitor data will be our outcome variable. We want to determine if we can predict air pollution levels based on other types of data, like road density and population density to see if we can use these data to predict air pollution in areas where there are no monitors.

Our outcome variable

The monitor data that we will be using comes from gravimetric monitors operated by the US Enivornmental Protection Agency (EPA). These monitors use a filtration system to specifically capture fine particulate matter. The weight of this matter is manually measured daily or weekly. See here for the EPA standard operating procedure for PM gravimetric analysis in 2008.

[source]

Here is an image of what the gravimetric monitors look like:

Gravimetric analysis is also used for emission testing. The same idea applies: a fresh filter is applied and the desired amount of time passes, then the filter is removed and weighed.

There are other monitoring systems that can provide hourly measurements, but we will not be using data from these monitors in our analysis. Gravimetric analysis is considered to be among the most accurate methods.

In our csv, the value column indicates the PM_2.5 monitor average for 2008 in mass of fine particles/volume of air for 876 gravimetric monitors. The units are micro grams of fine particulate mater (PM) that is less than 2.5 micrometers in diameter per cubic meter of air - mass concentration (ug/m³). Recall the the WHO exposure guideline is < 10 ug/m³ on average annually for PM_2.5.

Our predictor variables

There are 48 predictor variables with values for each of the 876 monitors included in our outcome variable. The data comes from the US Enivornmental Protection Agency (EPA), the National Aeronautics and Space Administration (NASA), the US Census, and the National Center for Health Statistics (NCHS).

Variable	Details
id	Monitor number – the county number is indicated before the decimal – the monitor number is indicated after the decimal Example: 1073.0023 is Jefferson county (1073) and .0023 one of 8 monitors
fips	Federal information processing standard number for the county where the monitor is located – 5 digit id code for counties (zero is often the first value and sometimes is not shown) – the first 2 numbers indicate the state – the last three numbers indicate the county Example: Alabama’s state code is 01 because it is first alphabetically (note: Alaska and Hawaii are not included because they are not part of the contiguous US)
Lat	Latitude of the monitor in degrees
Lon	Longitude of the monitor in degrees
state	State where the monitor is located
county	County where the monitor is located
city	City where the monitor is located
CMAQ	Estimated values of air pollution from a computational model called Community Multiscale Air Quality (CMAQ) – A monitoring system that simulates the physics of the atmosphere using chemistry and weather data to predict the air pollution – Does not use any of the PM2.5 gravimetric monitoring data. (There is a version that does use the gravimetric monitoring data, but not this one!) – Data from the EPA
zcta	Zip Code Tabulation Area where the monitor is located – Postal Zip codes are converted into “generalized areal representations” that are non-overlapping – Data from the 2010 Census
zcta_area	Land area of the zip code area in meters squared – Data from the 2010 Census
zcta_pop	Population in the zip code area – Data from the 2010 Census
imp_a500	Impervious surface measure – Within a circle with a radius of 500 meters around the monitor – Impervious surface are roads, concrete, parking lots, buildings – This is a measure of development
imp_a1000	Impervious surface measure – Within a circle with a radius of 1000 meters around the monitor
imp_a5000	Impervious surface measure – Within a circle with a radius of 5000 meters around the monitor
imp_a10000	Impervious surface measure – Within a circle with a radius of 10000 meters around the monitor
imp_a15000	Impervious surface measure – Within a circle with a radius of 15000 meters around the monitor
county_area	Land area of the county of the monitor in meters squared
county_pop	Population of the county of the monitor
Log_dist_to_prisec	Log (Natural log) distance to a primary or secondary road from the monitor – Highway or major road
log_pri_length_5000	Count of primary road length in meters in a circle with a radius of 5000 meters around the monitor (Natural log) – Highways only
log_pri_length_10000	Count of primary road length in meters in a circle with a radius of 10000 meters around the monitor (Natural log) – Highways only
log_pri_length_15000	Count of primary road length in meters in a circle with a radius of 15000 meters around the monitor (Natural log) – Highways only
log_pri_length_25000	Count of primary road length in meters in a circle with a radius of 25000 meters around the monitor (Natural log) – Highways only
log_prisec_length_500	Count of primary and secondary road length in meters in a circle with a radius of 500 meters around the monitor (Natural log) – Highway and secondary roads
log_prisec_length_1000	Count of primary and secondary road length in meters in a circle with a radius of 1000 meters around the monitor (Natural log) – Highway and secondary roads
log_prisec_length_5000	Count of primary and secondary road length in meters in a circle with a radius of 5000 meters around the monitor (Natural log) – Highway and secondary roads
log_prisec_length_10000	Count of primary and secondary road length in meters in a circle with a radius of 10000 meters around the monitor (Natural log) – Highway and secondary roads
log_prisec_length_15000	Count of primary and secondary road length in meters in a circle with a radius of 15000 meters around the monitor (Natural log) – Highway and secondary roads
log_prisec_length_25000	Count of primary and secondary road length in meters in a circle with a radius of 25000 meters around the monitor (Natural log) – Highway and secondary roads
log_nei_2008_pm25_sum_10000	Tons of emissions from major sources data base (annual data) sum of all sources within a circle with a radius of 10000 meters of distance around the monitor (Natural log)
log_nei_2008_pm25_sum_15000	Tons of emissions from major sources data base (annual data) sum of all sources within a circle with a radius of 15000 meters of distance around the monitor (Natural log)
log_nei_2008_pm25_sum_25000	Tons of emissions from major sources data base (annual data) sum of all sources within a circle with a radius of 25000 meters of distance around the monitor (Natural log)
log_nei_2008_pm10_sum_10000	Tons of emissions from major sources data base (annual data) sum of all sources within a circle with a radius of 10000 meters of distance around the monitor (Natural log)
log_nei_2008_pm10_sum_15000	Tons of emissions from major sources data base (annual data) sum of all sources within a circle with a radius of 15000 meters of distance around the monitor (Natural log)
log_nei_2008_pm10_sum_25000	Tons of emissions from major sources data base (annual data) sum of all sources within a circle with a radius of 25000 meters of distance around the monitor (Natural log)
popdens_county	Population density (number of people per kilometer squared area of the county)
popdens_zcta	Population density (number of people per kilometer squared area of zcta)
nohs	Percentage of people in zcta area where the monitor is that do not have a high school degree – Data from the Census
somehs	Percentage of people in zcta area where the monitor whose highest formal educational attainment was some high school education – Data from the Census
hs	Percentage of people in zcta area where the monitor whose highest formal educational attainment was completing a high school degree – Data from the Census
somecollege	Percentage of people in zcta area where the monitor whose highest formal educational attainment was completing some college education – Data from the Census
associate	Percentage of people in zcta area where the monitor whose highest formal educational attainment was completing an associate degree – Data from the Census
bachelor	Percentage of people in zcta area where the monitor whose highest formal educational attainment was a bachelor’s degree – Data from the Census
grad	Percentage of people in zcta area where the monitor whose highest formal educational attainment was a graduate degree – Data from the Census
pov	Percentage of people in zcta area where the monitor is that lived in poverty in 2008 - or would it have been 2007 guidelines??https://aspe.hhs.gov/2007-hhs-poverty-guidelines – Data from the Census
hs_orless	Percentage of people in zcta area where the monitor whose highest formal educational attainment was a high school degree or less (sum of nohs, somehs, and hs)
urc2013	2013 Urban-rural classification of the county where the monitor is located – 6 category variable - 1 is totally urban 6 is completely rural – Data from the National Center for Health Statistics](https://www.cdc.gov/nchs/index.htm){target="_blank"}
urc2006	2006 Urban-rural classification of the county where the monitor is located – 6 category variable - 1 is totally urban 6 is completely rural – Data from the National Center for Health Statistics
aod	Aerosol Optical Depth measurement from a NASA satellite – based on the defraction of a laser – used as a proxy of particulate pollution – unit-less - higher value indicates more pollution – Data from NASA

Many of these predictor variables have to do with the circular area around the monitor called the “buffer”. These are illustrated in the following figure:

[source]

Data Import

We have one CSV file that contains both our single outcome variable and all of our predictor variables.

Let’s import our data into R now so that we can explore the data further. We will call our data object pm for particulate matter.

pm <-readr::read_csv(here("docs", "pm25_data.csv"))

Data Exploration and Wrangling

The first step in performing a machine learning analysis is to explore the data to better understand the variables included in the data, as we may learn about important details about the data that we should keep in mind as we try to predict our outcome variable.

First let’s just get a general sense of our data. We can do that using the glimpse() function of the dplyr package (it is also in the tibble package).

We will also use the %>% pipe which can be used to define the input for later sequential steps. This will make more sense when we have multiple sequential steps using the same data object. To use the pipe notation we need to install and load dplyr as well.

For example here we will first grab the pm data object, then we use the glimpse() function on it based on the pipe notation.

pm %>%
  dplyr::glimpse()

Rows: 876
Columns: 50
$ id                          <dbl> 1003.001, 1027.000, 1033.100, 1049.100, 1…
$ value                       <dbl> 9.597647, 10.800000, 11.212174, 11.659091…
$ fips                        <dbl> 1003, 1027, 1033, 1049, 1055, 1069, 1073,…
$ lat                         <dbl> 30.49800, 33.28126, 34.75878, 34.28763, 3…
$ lon                         <dbl> -87.88141, -85.80218, -87.65056, -85.9683…
$ state                       <chr> "Alabama", "Alabama", "Alabama", "Alabama…
$ county                      <chr> "Baldwin", "Clay", "Colbert", "DeKalb", "…
$ city                        <chr> "Fairhope", "Ashland", "Muscle Shoals", "…
$ CMAQ                        <dbl> 8.098836, 9.766208, 9.402679, 8.534772, 9…
$ zcta                        <dbl> 36532, 36251, 35660, 35962, 35901, 36303,…
$ zcta_area                   <dbl> 190980522, 374132430, 16716984, 203836235…
$ zcta_pop                    <dbl> 27829, 5103, 9042, 8300, 20045, 30217, 90…
$ imp_a500                    <dbl> 0.01730104, 1.96972318, 19.17301038, 5.78…
$ imp_a1000                   <dbl> 1.4096021, 0.8531574, 11.1448962, 3.86764…
$ imp_a5000                   <dbl> 3.3360118, 0.9851479, 15.1786154, 1.23114…
$ imp_a10000                  <dbl> 1.9879187, 0.5208189, 9.7253870, 1.031646…
$ imp_a15000                  <dbl> 1.4386207, 0.3359198, 5.2472094, 0.973044…
$ county_area                 <dbl> 4117521611, 1564252280, 1534877333, 20126…
$ county_pop                  <dbl> 182265, 13932, 54428, 71109, 104430, 1015…
$ log_dist_to_prisec          <dbl> 4.648181, 7.219907, 5.760131, 3.721489, 5…
$ log_pri_length_5000         <dbl> 8.517193, 8.517193, 8.517193, 8.517193, 9…
$ log_pri_length_10000        <dbl> 9.210340, 9.210340, 9.274303, 10.409411, …
$ log_pri_length_15000        <dbl> 9.630228, 9.615805, 9.658899, 11.173626, …
$ log_pri_length_25000        <dbl> 11.32735, 10.12663, 10.15769, 11.90959, 1…
$ log_prisec_length_500       <dbl> 7.295356, 6.214608, 8.611945, 7.310155, 8…
$ log_prisec_length_1000      <dbl> 8.195119, 7.600902, 9.735569, 8.585843, 9…
$ log_prisec_length_5000      <dbl> 10.815042, 10.170878, 11.770407, 10.21420…
$ log_prisec_length_10000     <dbl> 11.88680, 11.40554, 12.84066, 11.50894, 1…
$ log_prisec_length_15000     <dbl> 12.205723, 12.042963, 13.282656, 12.35366…
$ log_prisec_length_25000     <dbl> 13.41395, 12.79980, 13.79973, 13.55979, 1…
$ log_nei_2008_pm25_sum_10000 <dbl> 0.318035438, 3.218632928, 6.573127301, 0.…
$ log_nei_2008_pm25_sum_15000 <dbl> 1.967358961, 3.218632928, 6.581917457, 3.…
$ log_nei_2008_pm25_sum_25000 <dbl> 5.067308, 3.218633, 6.875900, 4.887665, 4…
$ log_nei_2008_pm10_sum_10000 <dbl> 1.35588511, 3.31111648, 6.69187313, 0.000…
$ log_nei_2008_pm10_sum_15000 <dbl> 2.26783411, 3.31111648, 6.70127741, 3.350…
$ log_nei_2008_pm10_sum_25000 <dbl> 5.628728, 3.311116, 7.148858, 5.171920, 4…
$ popdens_county              <dbl> 44.265706, 8.906492, 35.460814, 35.330814…
$ popdens_zcta                <dbl> 145.716431, 13.639555, 540.887040, 40.718…
$ nohs                        <dbl> 3.3, 11.6, 7.3, 14.3, 4.3, 5.8, 7.1, 2.7,…
$ somehs                      <dbl> 4.9, 19.1, 15.8, 16.7, 13.3, 11.6, 17.1, …
$ hs                          <dbl> 25.1, 33.9, 30.6, 35.0, 27.8, 29.8, 37.2,…
$ somecollege                 <dbl> 19.7, 18.8, 20.9, 14.9, 29.2, 21.4, 23.5,…
$ associate                   <dbl> 8.2, 8.0, 7.6, 5.5, 10.1, 7.9, 7.3, 8.0, …
$ bachelor                    <dbl> 25.3, 5.5, 12.7, 7.9, 10.0, 13.7, 5.9, 17…
$ grad                        <dbl> 13.5, 3.1, 5.1, 5.8, 5.4, 9.8, 2.0, 8.7, …
$ pov                         <dbl> 6.1, 19.5, 19.0, 13.8, 8.8, 15.6, 25.5, 7…
$ hs_orless                   <dbl> 33.3, 64.6, 53.7, 66.0, 45.4, 47.2, 61.4,…
$ urc2013                     <dbl> 4, 6, 4, 6, 4, 4, 1, 1, 1, 1, 1, 1, 1, 2,…
$ urc2006                     <dbl> 5, 6, 4, 5, 4, 4, 1, 1, 1, 1, 1, 1, 1, 2,…
$ aod                         <dbl> 37.36364, 34.81818, 36.00000, 33.08333, 4…

We can see that there are 876 monitors and that we have 50 total variables - one of which is the outcome. In this case our outcome variable is called value.

Notice that some of the variables that we would think of as factors (categorical) are currently of class double as indicated by the <dbl> just to the right of the column names/variable names in the glimpse() output. For example the monitor ID (id), the Federal Information Processing Standard number for the county where the monitor was located (fips), as well as the zcta

Let’s convert these variables into factors. We can do this using the mutate_at() function of the dplyr package and the as.factor() base function.

In this case we are also using the magrittr assignment pipe or double pipe that looks like this %<>% of the magrittr package. This allows us use the pm data as input but also reassign the output to the same data object name.

pm %<>%
  dplyr::mutate_at(vars(id, fips, zcta), as.factor) 

glimpse(pm)

Rows: 876
Columns: 50
$ id                          <fct> 1003.001, 1027.0001, 1033.1002, 1049.1003…
$ value                       <dbl> 9.597647, 10.800000, 11.212174, 11.659091…
$ fips                        <fct> 1003, 1027, 1033, 1049, 1055, 1069, 1073,…
$ lat                         <dbl> 30.49800, 33.28126, 34.75878, 34.28763, 3…
$ lon                         <dbl> -87.88141, -85.80218, -87.65056, -85.9683…
$ state                       <chr> "Alabama", "Alabama", "Alabama", "Alabama…
$ county                      <chr> "Baldwin", "Clay", "Colbert", "DeKalb", "…
$ city                        <chr> "Fairhope", "Ashland", "Muscle Shoals", "…
$ CMAQ                        <dbl> 8.098836, 9.766208, 9.402679, 8.534772, 9…
$ zcta                        <fct> 36532, 36251, 35660, 35962, 35901, 36303,…
$ zcta_area                   <dbl> 190980522, 374132430, 16716984, 203836235…
$ zcta_pop                    <dbl> 27829, 5103, 9042, 8300, 20045, 30217, 90…
$ imp_a500                    <dbl> 0.01730104, 1.96972318, 19.17301038, 5.78…
$ imp_a1000                   <dbl> 1.4096021, 0.8531574, 11.1448962, 3.86764…
$ imp_a5000                   <dbl> 3.3360118, 0.9851479, 15.1786154, 1.23114…
$ imp_a10000                  <dbl> 1.9879187, 0.5208189, 9.7253870, 1.031646…
$ imp_a15000                  <dbl> 1.4386207, 0.3359198, 5.2472094, 0.973044…
$ county_area                 <dbl> 4117521611, 1564252280, 1534877333, 20126…
$ county_pop                  <dbl> 182265, 13932, 54428, 71109, 104430, 1015…
$ log_dist_to_prisec          <dbl> 4.648181, 7.219907, 5.760131, 3.721489, 5…
$ log_pri_length_5000         <dbl> 8.517193, 8.517193, 8.517193, 8.517193, 9…
$ log_pri_length_10000        <dbl> 9.210340, 9.210340, 9.274303, 10.409411, …
$ log_pri_length_15000        <dbl> 9.630228, 9.615805, 9.658899, 11.173626, …
$ log_pri_length_25000        <dbl> 11.32735, 10.12663, 10.15769, 11.90959, 1…
$ log_prisec_length_500       <dbl> 7.295356, 6.214608, 8.611945, 7.310155, 8…
$ log_prisec_length_1000      <dbl> 8.195119, 7.600902, 9.735569, 8.585843, 9…
$ log_prisec_length_5000      <dbl> 10.815042, 10.170878, 11.770407, 10.21420…
$ log_prisec_length_10000     <dbl> 11.88680, 11.40554, 12.84066, 11.50894, 1…
$ log_prisec_length_15000     <dbl> 12.205723, 12.042963, 13.282656, 12.35366…
$ log_prisec_length_25000     <dbl> 13.41395, 12.79980, 13.79973, 13.55979, 1…
$ log_nei_2008_pm25_sum_10000 <dbl> 0.318035438, 3.218632928, 6.573127301, 0.…
$ log_nei_2008_pm25_sum_15000 <dbl> 1.967358961, 3.218632928, 6.581917457, 3.…
$ log_nei_2008_pm25_sum_25000 <dbl> 5.067308, 3.218633, 6.875900, 4.887665, 4…
$ log_nei_2008_pm10_sum_10000 <dbl> 1.35588511, 3.31111648, 6.69187313, 0.000…
$ log_nei_2008_pm10_sum_15000 <dbl> 2.26783411, 3.31111648, 6.70127741, 3.350…
$ log_nei_2008_pm10_sum_25000 <dbl> 5.628728, 3.311116, 7.148858, 5.171920, 4…
$ popdens_county              <dbl> 44.265706, 8.906492, 35.460814, 35.330814…
$ popdens_zcta                <dbl> 145.716431, 13.639555, 540.887040, 40.718…
$ nohs                        <dbl> 3.3, 11.6, 7.3, 14.3, 4.3, 5.8, 7.1, 2.7,…
$ somehs                      <dbl> 4.9, 19.1, 15.8, 16.7, 13.3, 11.6, 17.1, …
$ hs                          <dbl> 25.1, 33.9, 30.6, 35.0, 27.8, 29.8, 37.2,…
$ somecollege                 <dbl> 19.7, 18.8, 20.9, 14.9, 29.2, 21.4, 23.5,…
$ associate                   <dbl> 8.2, 8.0, 7.6, 5.5, 10.1, 7.9, 7.3, 8.0, …
$ bachelor                    <dbl> 25.3, 5.5, 12.7, 7.9, 10.0, 13.7, 5.9, 17…
$ grad                        <dbl> 13.5, 3.1, 5.1, 5.8, 5.4, 9.8, 2.0, 8.7, …
$ pov                         <dbl> 6.1, 19.5, 19.0, 13.8, 8.8, 15.6, 25.5, 7…
$ hs_orless                   <dbl> 33.3, 64.6, 53.7, 66.0, 45.4, 47.2, 61.4,…
$ urc2013                     <dbl> 4, 6, 4, 6, 4, 4, 1, 1, 1, 1, 1, 1, 1, 2,…
$ urc2006                     <dbl> 5, 6, 4, 5, 4, 4, 1, 1, 1, 1, 1, 1, 1, 2,…
$ aod                         <dbl> 37.36364, 34.81818, 36.00000, 33.08333, 4…

Great! Now we can see that these variables are now factors as indicated by <fct> after the variable name.

Packages to get a sense of the data

The skim() function of the skimr package is also really helpful for getting a general sense of your data.

skim(pm)

Data summary

Name	pm
Number of rows	876
Number of columns	50
_______________________
Column type frequency:
character	3
factor	3
numeric	44
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
state	0	1	4	20	0	49	0
county	0	1	3	20	0	471	0
city	0	1	4	48	0	607	0

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
id	0	1	FALSE	876	100: 1, 102: 1, 103: 1, 104: 1
fips	0	1	FALSE	569	170: 12, 603: 10, 261: 9, 107: 8
zcta	0	1	FALSE	842	475: 3, 110: 2, 160: 2, 290: 2

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
value	0	1	10.81	2.580000e+00	3.02	9.27	11.15	12.37	2.316000e+01	▂▆▇▁▁
lat	0	1	38.48	4.620000e+00	25.47	35.03	39.30	41.66	4.840000e+01	▁▃▅▇▂
lon	0	1	-91.74	1.496000e+01	-124.18	-99.16	-87.47	-80.69	-6.804000e+01	▃▂▃▇▃
CMAQ	0	1	8.41	2.970000e+00	1.63	6.53	8.62	10.24	2.313000e+01	▃▇▃▁▁
zcta_area	0	1	183173481.91	5.425989e+08	15459.00	14204601.75	37653560.50	160041508.25	8.164821e+09	▇▁▁▁▁
zcta_pop	0	1	24227.58	1.777216e+04	0.00	9797.00	22014.00	35004.75	9.539700e+04	▇▇▃▁▁
imp_a500	0	1	24.72	1.934000e+01	0.00	3.70	25.12	40.22	6.961000e+01	▇▅▆▃▂
imp_a1000	0	1	24.26	1.802000e+01	0.00	5.32	24.53	38.59	6.750000e+01	▇▅▆▃▁
imp_a5000	0	1	19.93	1.472000e+01	0.05	6.79	19.07	30.11	7.460000e+01	▇▆▃▁▁
imp_a10000	0	1	15.82	1.381000e+01	0.09	4.54	12.36	24.17	7.209000e+01	▇▃▂▁▁
imp_a15000	0	1	13.43	1.312000e+01	0.11	3.24	9.67	20.55	7.110000e+01	▇▃▁▁▁
county_area	0	1	3768701992.12	6.212830e+09	33703512.00	1116536297.50	1690826566.50	2878192209.00	5.194723e+10	▇▁▁▁▁
county_pop	0	1	687298.44	1.293489e+06	783.00	100948.00	280730.50	743159.00	9.818605e+06	▇▁▁▁▁
log_dist_to_prisec	0	1	6.19	1.410000e+00	-1.46	5.43	6.36	7.15	1.045000e+01	▁▁▃▇▁
log_pri_length_5000	0	1	9.82	1.080000e+00	8.52	8.52	10.05	10.73	1.205000e+01	▇▂▆▅▂
log_pri_length_10000	0	1	10.92	1.130000e+00	9.21	9.80	11.17	11.83	1.302000e+01	▇▂▇▇▃
log_pri_length_15000	0	1	11.50	1.150000e+00	9.62	10.87	11.72	12.40	1.359000e+01	▆▂▇▇▃
log_pri_length_25000	0	1	12.24	1.100000e+00	10.13	11.69	12.46	13.05	1.436000e+01	▅▃▇▇▃
log_prisec_length_500	0	1	6.99	9.500000e-01	6.21	6.21	6.21	7.82	9.400000e+00	▇▁▂▂▁
log_prisec_length_1000	0	1	8.56	7.900000e-01	7.60	7.60	8.66	9.20	1.047000e+01	▇▅▆▃▁
log_prisec_length_5000	0	1	11.28	7.800000e-01	8.52	10.91	11.42	11.83	1.278000e+01	▁▁▃▇▃
log_prisec_length_10000	0	1	12.41	7.300000e-01	9.21	11.99	12.53	12.94	1.385000e+01	▁▁▃▇▅
log_prisec_length_15000	0	1	13.03	7.200000e-01	9.62	12.59	13.13	13.57	1.441000e+01	▁▁▃▇▅
log_prisec_length_25000	0	1	13.82	7.000000e-01	10.13	13.38	13.92	14.35	1.523000e+01	▁▁▃▇▆
log_nei_2008_pm25_sum_10000	0	1	3.97	2.350000e+00	0.00	2.15	4.29	5.69	9.120000e+00	▆▅▇▆▂
log_nei_2008_pm25_sum_15000	0	1	4.72	2.250000e+00	0.00	3.47	5.00	6.35	9.420000e+00	▃▃▇▇▂
log_nei_2008_pm25_sum_25000	0	1	5.67	2.110000e+00	0.00	4.66	5.91	7.28	9.650000e+00	▂▂▇▇▃
log_nei_2008_pm10_sum_10000	0	1	4.35	2.320000e+00	0.00	2.69	4.62	6.07	9.340000e+00	▅▅▇▇▂
log_nei_2008_pm10_sum_15000	0	1	5.10	2.180000e+00	0.00	3.87	5.39	6.72	9.710000e+00	▂▃▇▇▂
log_nei_2008_pm10_sum_25000	0	1	6.07	2.010000e+00	0.00	5.10	6.37	7.52	9.880000e+00	▁▂▆▇▃
popdens_county	0	1	551.76	1.711510e+03	0.26	40.77	156.67	510.81	2.682191e+04	▇▁▁▁▁
popdens_zcta	0	1	1279.66	2.757490e+03	0.00	101.15	610.35	1382.52	3.041884e+04	▇▁▁▁▁
nohs	0	1	6.99	7.210000e+00	0.00	2.70	5.10	8.80	1.000000e+02	▇▁▁▁▁
somehs	0	1	10.17	6.200000e+00	0.00	5.90	9.40	13.90	7.220000e+01	▇▂▁▁▁
hs	0	1	30.32	1.140000e+01	0.00	23.80	30.75	36.10	1.000000e+02	▂▇▂▁▁
somecollege	0	1	21.58	8.600000e+00	0.00	17.50	21.30	24.70	1.000000e+02	▆▇▁▁▁
associate	0	1	7.13	4.010000e+00	0.00	4.90	7.10	8.80	7.140000e+01	▇▁▁▁▁
bachelor	0	1	14.90	9.710000e+00	0.00	8.80	12.95	19.22	1.000000e+02	▇▂▁▁▁
grad	0	1	8.91	8.650000e+00	0.00	3.90	6.70	11.00	1.000000e+02	▇▁▁▁▁
pov	0	1	14.95	1.133000e+01	0.00	6.50	12.10	21.22	6.590000e+01	▇▅▂▁▁
hs_orless	0	1	47.48	1.675000e+01	0.00	37.92	48.65	59.10	1.000000e+02	▁▃▇▃▁
urc2013	0	1	2.92	1.520000e+00	1.00	2.00	3.00	4.00	6.000000e+00	▇▅▃▂▁
urc2006	0	1	2.97	1.520000e+00	1.00	2.00	3.00	4.00	6.000000e+00	▇▅▃▂▁
aod	0	1	43.70	1.956000e+01	5.00	31.66	40.17	49.67	1.430000e+02	▃▇▁▁▁

Notice how there is a column called n_missing about the number of values that are missing. It looks like our data is very complete and we do not have any missing data. This is also indicated by the complete_rate variable, which shows the ratio of completeness, in our case all variables have a value of 1 indicating they are fully complete.

The n_unqiue column shows us the number of unique values for each of our columns. We can see that there are 49 states represented in the data, and we know that the data should be of the contiguous states. Let’s take a look to see which states are included:

# Scroll through the output!
pm %>% 
  distinct(state) %>%
  print(n = 1e3)

# A tibble: 49 x 1
   state               
   <chr>               
 1 Alabama             
 2 Arizona             
 3 Arkansas            
 4 California          
 5 Colorado            
 6 Connecticut         
 7 Delaware            
 8 District Of Columbia
 9 Florida             
10 Georgia             
11 Idaho               
12 Illinois            
13 Indiana             
14 Iowa                
15 Kansas              
16 Kentucky            
17 Louisiana           
18 Maine               
19 Maryland            
20 Massachusetts       
21 Michigan            
22 Minnesota           
23 Mississippi         
24 Missouri            
25 Montana             
26 Nebraska            
27 Nevada              
28 New Hampshire       
29 New Jersey          
30 New Mexico          
31 New York            
32 North Carolina      
33 North Dakota        
34 Ohio                
35 Oklahoma            
36 Oregon              
37 Pennsylvania        
38 Rhode Island        
39 South Carolina      
40 South Dakota        
41 Tennessee           
42 Texas               
43 Utah                
44 Vermont             
45 Virginia            
46 Washington          
47 West Virginia       
48 Wisconsin           
49 Wyoming             

Looks like “District of Columbia” is being included as a state. We can see that indeed Alaska and Hawaii are not included in the data.

Here is another method of looking at the data using the dfSummary() function of the summarytoolspackage. We need to copy and paste the output into the rmarkdown.

dfSummary(pm, plain.ascii = FALSE, style = "grid", 
          graph.magnif = 0.45,  tmp.img.dir = "tmp")

Dimensions: 876 x 50
Duplicates: 0

No	Variable	Stats / Values	Freqs (% of Valid)	Graph	Valid	Missing
1	id [factor]	1. 1003.001 2. 1027.0001 3. 1033.1002 4. 1049.1003 5. 1055.001 6. 1069.0003 7. 1073.0023 8. 1073.1005 9. 1073.1009 10. 1073.101 [ 866 others ]	1 ( 0.1%) 1 ( 0.1%) 1 ( 0.1%) 1 ( 0.1%) 1 ( 0.1%) 1 ( 0.1%) 1 ( 0.1%) 1 ( 0.1%) 1 ( 0.1%) 1 ( 0.1%) 866 (98.9%)		876 (100%)	0 (0%)
2	value [numeric]	Mean (sd) : 10.8 (2.6) min < med < max: 3 < 11.2 < 23.2 IQR (CV) : 3.1 (0.2)	875 distinct values		876 (100%)	0 (0%)
3	fips [factor]	1. 1003 2. 1027 3. 1033 4. 1049 5. 1055 6. 1069 7. 1073 8. 1089 9. 1097 10. 1101 [ 559 others ]	1 ( 0.1%) 1 ( 0.1%) 1 ( 0.1%) 1 ( 0.1%) 1 ( 0.1%) 1 ( 0.1%) 8 ( 0.9%) 1 ( 0.1%) 2 ( 0.2%) 1 ( 0.1%) 858 (98.0%)		876 (100%)	0 (0%)
4	lat [numeric]	Mean (sd) : 38.5 (4.6) min < med < max: 25.5 < 39.3 < 48.4 IQR (CV) : 6.6 (0.1)	876 distinct values		876 (100%)	0 (0%)
5	lon [numeric]	Mean (sd) : -91.7 (15) min < med < max: -124.2 < -87.5 < -68 IQR (CV) : 18.5 (-0.2)	876 distinct values		876 (100%)	0 (0%)
6	state [character]	1. California 2. Ohio 3. Illinois 4. Indiana 5. North Carolina 6. Pennsylvania 7. Michigan 8. Florida 9. Georgia 10. Texas [ 39 others ]	85 ( 9.7%) 44 ( 5.0%) 38 ( 4.3%) 36 ( 4.1%) 35 ( 4.0%) 32 ( 3.7%) 30 ( 3.4%) 29 ( 3.3%) 28 ( 3.2%) 27 ( 3.1%) 492 (56.2%)		876 (100%)	0 (0%)
7	county [character]	1. Jefferson 2. Cook 3. Hamilton 4. Lake 5. Los Angeles 6. Wayne 7. Washington 8. Cuyahoga 9. Jackson 10. Madison [ 461 others ]	18 ( 2.1%) 12 ( 1.4%) 11 ( 1.3%) 11 ( 1.3%) 10 ( 1.1%) 10 ( 1.1%) 9 ( 1.0%) 7 ( 0.8%) 7 ( 0.8%) 7 ( 0.8%) 774 (88.4%)		876 (100%)	0 (0%)
8	city [character]	1. Not in a city 2. New York 3. Cleveland 4. Baltimore 5. Chicago 6. Detroit 7. Milwaukee 8. New Haven 9. Philadelphia 10. Springfield [ 597 others ]	103 (11.8%) 9 ( 1.0%) 6 ( 0.7%) 5 ( 0.6%) 5 ( 0.6%) 5 ( 0.6%) 5 ( 0.6%) 5 ( 0.6%) 5 ( 0.6%) 5 ( 0.6%) 723 (82.5%)		876 (100%)	0 (0%)
9	CMAQ  | Me [numeric]	an (sd) : 8.4 (3)  | 60 min < med < max: 1.6 < 8.6 < 23.1 IQR (CV) : 3.7 (0.4)	1 distinct values |  | 87	6  | 0 (100%)	(0%)
10	zcta [factor]	1. 1022 2. 1103 3. 1201 4. 1608 5. 1832 6. 1840 7. 1863 8. 1904 9. 2113 10. 2119 [ 832 others ]	1 ( 0.1%) 2 ( 0.2%) 1 ( 0.1%) 2 ( 0.2%) 1 ( 0.1%) 1 ( 0.1%) 1 ( 0.1%) 1 ( 0.1%) 1 ( 0.1%) 1 ( 0.1%) 864 (98.6%)		876 (100%)	0 (0%)
11	zcta_area [numeric]	Mean (sd) : 183173481.9 (542598878.5) min < med < max: 15459 < 37653560.5 < 8164820625 IQR (CV) : 145836906.5 (3)	842 distinct values		876 (100%)	0 (0%)
12	zcta_pop [numeric]	Mean (sd) : 24227.6 (17772.2) min < med < max: 0 < 22014 < 95397 IQR (CV) : 25207.8 (0.7)	837 distinct values		876 (100%)	0 (0%)
13	imp_a500 [numeric]	Mean (sd) : 24.7 (19.3) min < med < max: 0 < 25.1 < 69.6 IQR (CV) : 36.5 (0.8)	816 distinct values		876 (100%)	0 (0%)
14	imp_a1000 [numeric]	Mean (sd) : 24.3 (18) min < med < max: 0 < 24.5 < 67.5 IQR (CV) : 33.3 (0.7)	860 distinct values		876 (100%)	0 (0%)
15	imp_a5000 [numeric]	Mean (sd) : 19.9 (14.7) min < med < max: 0.1 < 19.1 < 74.6 IQR (CV) : 23.3 (0.7)	870 distinct values		876 (100%)	0 (0%)
16	imp_a10000 [numeric]	Mean (sd) : 15.8 (13.8) min < med < max: 0.1 < 12.4 < 72.1 IQR (CV) : 19.6 (0.9)	870 distinct values		876 (100%)	0 (0%)
17	imp_a15000 [numeric]	Mean (sd) : 13.4 (13.1) min < med < max: 0.1 < 9.7 < 71.1 IQR (CV) : 17.3 (1)	870 distinct values		876 (100%)	0 (0%)
18	county_area [numeric]	Mean (sd) : 3768701992.1 (6212829553.6) min < med < max: 33703512 < 1690826566.5 < 51947229509 IQR (CV) : 1761655911.5 (1.6)	564 distinct values		876 (100%)	0 (0%)
19	county_pop [numeric]	Mean (sd) : 687298.4 (1293488.7) min < med < max: 783 < 280730.5 < 9818605 IQR (CV) : 642211 (1.9)	564 distinct values		876 (100%)	0 (0%)
20	log_dist_to_prisec [numeric]	Mean (sd) : 6.2 (1.4) min < med < max: -1.5 < 6.4 < 10.5 IQR (CV) : 1.7 (0.2)	870 distinct values		876 (100%)	0 (0%)
21	log_pri_length_5000 [numeric]	Mean (sd) : 9.8 (1.1) min < med < max: 8.5 < 10.1 < 12 IQR (CV) : 2.2 (0.1)	586 distinct values		876 (100%)	0 (0%)
22	log_pri_length_10000 [numeric]	Mean (sd) : 10.9 (1.1) min < med < max: 9.2 < 11.2 < 13 IQR (CV) : 2 (0.1)	687 distinct values		876 (100%)	0 (0%)
23	log_pri_length_15000 [numeric]	Mean (sd) : 11.5 (1.1) min < med < max: 9.6 < 11.7 < 13.6 IQR (CV) : 1.5 (0.1)	726 distinct values		876 (100%)	0 (0%)
24	log_pri_length_25000 [numeric]	Mean (sd) : 12.2 (1.1) min < med < max: 10.1 < 12.5 < 14.4 IQR (CV) : 1.4 (0.1)	787 distinct values		876 (100%)	0 (0%)
25	log_prisec_length_500 [numeric]	Mean (sd) : 7 (1) min < med < max: 6.2 < 6.2 < 9.4 IQR (CV) : 1.6 (0.1)	382 distinct values		876 (100%)	0 (0%)
26	log_prisec_length_1000 [numeric]	Mean (sd) : 8.6 (0.8) min < med < max: 7.6 < 8.7 < 10.5 IQR (CV) : 1.6 (0.1)	591 distinct values		876 (100%)	0 (0%)
27	log_prisec_length_5000 [numeric]	Mean (sd) : 11.3 (0.8) min < med < max: 8.5 < 11.4 < 12.8 IQR (CV) : 0.9 (0.1)	852 distinct values		876 (100%)	0 (0%)
28	log_prisec_length_10000 [numeric]	Mean (sd) : 12.4 (0.7) min < med < max: 9.2 < 12.5 < 13.8 IQR (CV) : 1 (0.1)	867 distinct values		876 (100%)	0 (0%)
29	log_prisec_length_15000 [numeric]	Mean (sd) : 13 (0.7) min < med < max: 9.6 < 13.1 < 14.4 IQR (CV) : 1 (0.1)	869 distinct values		876 (100%)	0 (0%)
30	log_prisec_length_25000 [numeric]	Mean (sd) : 13.8 (0.7) min < med < max: 10.1 < 13.9 < 15.2 IQR (CV) : 1 (0.1)	870 distinct values		876 (100%)	0 (0%)
31	log_nei_2008_pm25_sum_10000 [numeric]	Mean (sd) : 4 (2.4) min < med < max: 0 < 4.3 < 9.1 IQR (CV) : 3.5 (0.6)	828 distinct values		876 (100%)	0 (0%)
32	log_nei_2008_pm25_sum_15000 [numeric]	Mean (sd) : 4.7 (2.2) min < med < max: 0 < 5 < 9.4 IQR (CV) : 2.9 (0.5)	855 distinct values		876 (100%)	0 (0%)
33	log_nei_2008_pm25_sum_25000 [numeric]	Mean (sd) : 5.7 (2.1) min < med < max: 0 < 5.9 < 9.7 IQR (CV) : 2.6 (0.4)	860 distinct values		876 (100%)	0 (0%)
34	log_nei_2008_pm10_sum_10000 [numeric]	Mean (sd) : 4.3 (2.3) min < med < max: 0 < 4.6 < 9.3 IQR (CV) : 3.4 (0.5)	829 distinct values		876 (100%)	0 (0%)
35	log_nei_2008_pm10_sum_15000 [numeric]	Mean (sd) : 5.1 (2.2) min < med < max: 0 < 5.4 < 9.7 IQR (CV) : 2.8 (0.4)	855 distinct values		876 (100%)	0 (0%)
36	log_nei_2008_pm10_sum_25000 [numeric]	Mean (sd) : 6.1 (2) min < med < max: 0 < 6.4 < 9.9 IQR (CV) : 2.4 (0.3)	860 distinct values		876 (100%)	0 (0%)
37	popdens_county [numeric]	Mean (sd) : 551.8 (1711.5) min < med < max: 0.3 < 156.7 < 26821.9 IQR (CV) : 470 (3.1)	564 distinct values		876 (100%)	0 (0%)
38	popdens_zcta [numeric]	Mean (sd) : 1279.7 (2757.5) min < med < max: 0 < 610.3 < 30418.8 IQR (CV) : 1281.4 (2.2)	840 distinct values		876 (100%)	0 (0%)
39	nohs [numeric]	Mean (sd) : 7 (7.2) min < med < max: 0 < 5.1 < 100 IQR (CV) : 6.1 (1)	215 distinct values		876 (100%)	0 (0%)
40	somehs [numeric]	Mean (sd) : 10.2 (6.2) min < med < max: 0 < 9.4 < 72.2 IQR (CV) : 8 (0.6)	230 distinct values		876 (100%)	0 (0%)
41	hs [numeric]	Mean (sd) : 30.3 (11.4) min < med < max: 0 < 30.8 < 100 IQR (CV) : 12.3 (0.4)	347 distinct values		876 (100%)	0 (0%)
42	somecollege [numeric]	Mean (sd) : 21.6 (8.6) min < med < max: 0 < 21.3 < 100 IQR (CV) : 7.2 (0.4)	240 distinct values		876 (100%)	0 (0%)
43	associate [numeric]	Mean (sd) : 7.1 (4) min < med < max: 0 < 7.1 < 71.4 IQR (CV) : 3.9 (0.6)	157 distinct values		876 (100%)	0 (0%)
44	bachelor [numeric]	Mean (sd) : 14.9 (9.7) min < med < max: 0 < 12.9 < 100 IQR (CV) : 10.4 (0.7)	301 distinct values		876 (100%)	0 (0%)
45	grad [numeric]	Mean (sd) : 8.9 (8.6) min < med < max: 0 < 6.7 < 100 IQR (CV) : 7.1 (1)	245 distinct values		876 (100%)	0 (0%)
46	pov [numeric]	Mean (sd) : 15 (11.3) min < med < max: 0 < 12.1 < 65.9 IQR (CV) : 14.7 (0.8)	345 distinct values		876 (100%)	0 (0%)
47	hs_orless [numeric]	Mean (sd) : 47.5 (16.8) min < med < max: 0 < 48.7 < 100 IQR (CV) : 21.2 (0.4)	464 distinct values		876 (100%)	0 (0%)
48	urc2013 [numeric]	Mean (sd) : 2.9 (1.5) min < med < max: 1 < 3 < 6 IQR (CV) : 2 (0.5)	1 : 203 (23.2%) 2 : 163 (18.6%) 3 : 228 (26.0%) 4 : 123 (14.0%) 5 : 101 (11.5%) 6 : 58 ( 6.6%)		876 (100%)	0 (0%)
49	urc2006 [numeric]	Mean (sd) : 3 (1.5) min < med < max: 1 < 3 < 6 IQR (CV) : 2 (0.5)	1 : 195 (22.3%) 2 : 162 (18.5%) 3 : 221 (25.2%) 4 : 127 (14.5%) 5 : 115 (13.1%) 6 : 56 ( 6.4%)		876 (100%)	0 (0%)
50	aod [numeric]	Mean (sd) : 43.7 (19.6) min < med < max: 5 < 40.2 < 143 IQR (CV) : 18 (0.4)	581 distinct values		876 (100%)	0 (0%)

We can see that for many variables there are many low values as the distribution shows two peaks, one near zero and another with a higher value. This is true for the imp variables (measures of development), the nei variables (measures of emission sources) and the road density variables. We can also see that the range of some of the variables is very large, in particular the area and population related variables.

Evaluate correlation among possible predictors

In prediction analyses, it is also useful to evaluate if any of the variables are correlated.

Intuitively we can expect some of our variables to be correlated.

Let’s first take a look at all of our numeric variabels with thecorrplot package: The corrplot package is another option to look at correlation among possible predictors. This is a great option if we have many predictors. First we need to create a correlation matrix using the cor() function of the stats package (which is loaded automatically).

#library(RColorBrewer) #need this package for the color
PM_cor <- cor(pm %>% dplyr::select_if(is.numeric))
corrplot::corrplot(PM_cor, tl.cex = 0.5)

corrplot(abs(PM_cor), order = "AOE", tl.cex = 0.5,  cl.lim = c(0, 1))

corrplot(PM_cor, diag = FALSE, order = "FPC",
         tl.pos = "td", tl.cex = 0.5, method = "color", type = "upper", col = brewer.pal(n = 8, name = "PuOr"))

corrplot(PM_cor, diag = FALSE,
         tl.pos = "td", tl.cex = 0.5, method = "color", type = "upper")

Using ggcorplot

library(ggcorrplot)
ggcorrplot(PM_cor, hc.order = TRUE, type = "lower", tl.cex = 5)

We can see that the the development variables (imp) variables are correlated with each other as we might expect. We also see that the road density variables seem to be correlated with each other, and the emission variables seem to be correlated with each other. We can take a closer look using the ggcorr() function and the ggpairs() function of the GGally package. To select our variables of interest we can use the select() function with the contains() function of the tidyr package.

First let’s look at the imp/development variables.

select(pm, contains("imp")) %>%
  ggcorr(palette = "RdBu", label = TRUE)

select(pm, contains("imp")) %>%
  ggpairs()

Indeed, we can see that imp_a1000 and imp_a500 are perfectly correlated, as well as imp_a10000, imp_a15000.

Now let’s take a look at the road density data:

select(pm, contains("pri")) %>%
  ggcorr(palette = "RdBu",  hjust = .85, size = 3,
       layout.exp=2, label = TRUE)

select(pm, contains("pri")) %>%
  ggpairs()

We can see that while some of the road density variables are highly correlated with one another (or anti-correlated), there arr some variables that are less correlated with one another.

Finally let’s look at the emission variables.

select(pm, contains("nei")) %>%
  ggcorr(palette = "RdBu",  hjust = .85, size = 3,
       layout.exp=2, label = TRUE)

select(pm, contains("nei")) %>%
  ggpairs()

We would also expect the population density data might correlate with some of these variables. Let’s take a look.

pm %>%
select(log_nei_2008_pm25_sum_10000, popdens_county, log_pri_length_10000, imp_a10000) %>%
  ggcorr(palette = "RdBu",  hjust = .85, size = 3,
       layout.exp=2, label = TRUE)

pm %>%
select(log_nei_2008_pm25_sum_10000, popdens_county, log_pri_length_10000, imp_a10000, county_pop) %>%
  ggpairs()

Interesting, so these variables don’t appear to be highly correlated, therefore we might need variables from each of the categories to predict our monitor PM_2.5 pollution values.

We seem to have some pretty extreme population values though, so let’s see what happens when we take the log value.

pm %>%
  mutate(log_popdens_county= log(popdens_county)) %>%
select(log_nei_2008_pm25_sum_10000, log_popdens_county, log_pri_length_10000, imp_a10000) %>%
  ggcorr(palette = "RdBu",  hjust = .85, size = 3,
       layout.exp=2, label = TRUE)

pm %>%
  mutate(log_popdens_county= log(popdens_county)) %>%
  mutate(log_pop_county = log(county_pop)) %>%
select(log_nei_2008_pm25_sum_10000, log_popdens_county, log_pri_length_10000, imp_a10000, log_pop_county) %>%
  ggpairs()

Indeed this increased the correlation, but variables from each of these categories may still prove to be useful for prediction.

Data Analysis

Now that we have a sense of what our data is like we can get started with data analysis.

The machine learning process

There are two major types of machine learning:

1. Unsupervised
   
2. Supervised

Unsupervised learning is used to learn about the structure of the data without knowing much about the data. We let the data reveal properties about itself. Examples include clustering the data into groups or reducing the dimensionality of the data using methods like principal component analysis (which we will describe in more detail later) to capture patterns of varaince within the data.

[source]

In contrast, in supervised learning we have some knowlege about the data that we want to use to create a model to be able to generalize about other similar data.

There are two distinct goals of supervised machine learning:

1. Prediction
   
2. Classification

[source]

We will be performing a prediction analysis (which is also refered to as regression), which aims to predict continuous outcome variables given a number of predictors/explanatory variables/features/parameters, as we have already described.

Classification on the other hand aims to discern or predict group identity for a categorical outcome based on a number of predictors/explanatory variables/features/parameters.

The overall process is the same in either case and involes the following steps (which will each be explained in detail):

1. Data exploration

We have already performed this step to get a sense of the data. It is important to know if we have NA values, to understand the class of variables, and if to determine if there are any redudant variables that might need to removed.

2. Data splitting

The data needs to be split into two peices: a training set and a testing set. The training set will be used to optimize the model, while the testing set will be used to evaluate model performance.

3. Variable assignment and pre-processing

Both the training and testing data needs to be processed so that the data is compatible and optimized to be used with the model. This involves assigning variables to specific roles within the model and preprocessing like scaling variables and removing redudant variables.

4. Model specification, fitting, tuning and performance evaluation using the training data

The model needs to first be fit to the training data. First the method or algorithm in which the model will be fit is specified (regression, random forests etc.). Then in both classification and prediction, the model is fit to the traing data and the explanatory variables are used to estimate numeric values (in the case of prediction) or categorical values (in the case of classification) of the outcome variable of interest. If the model fits well than these estimated values will be very similar to the true outcome variable values. If the model does not fit well, than these estimates will be more dissimal from the true outcome variable. In this case, aspects about the model may need to be modified to improve the similarity of the estimates with that ofthe true outcome values. One way to optimize model performance is a process called tuning in which different model hyperparameter options are tested to determine the best option for model performance.

5. Overall model performance evaluation

Model performance is assessed as the similarity between the estimates of the outcome variable produced by the model and the true outcome variable values. This is done typically as an interative process with the training data along side modification of the model untill the performance using the training data is satisfactory. At this point, the final model performance is assessed using the testing data. This then gives an estimate about how well the model will predict or classify the outcome variable of interest with new independent data. Ideally one would also perform an evaluation with independent data to provide a sense of how generalizable the model is to other data sources.

The tidymodels ecosystem

To perform our analysis we will be using the tidymodels suite of packages. You may be familiar with the older packages caret or mlr which are also for machine learning and modeling but are not a part of the tidyverse. Max Kuhn describes tidymodels like this:

“Other packages, such as caret and mlr, help to solve the R model API issue. These packages do a lot of other things too: preprocessing, model tuning, resampling, feature selection, ensembling, and so on. In the tidyverse, we strive to make our packages modular and parsnip is designed only to solve the interface issue. It is not designed to be a drop-in replacement for caret. The tidymodels package collection, which includes parsnip, has other packages for many of these tasks, and they are designed to work together. We are working towards higher-level APIs that can replicate and extend what the current model packages can do.”

There are many packages in the tidymodels ecosystem which assist with the various steps of the machine learning process:

This is a depiction of how these tools help perform the overall machine learning process:

The major benefits of tidymodels

1. Standardized workflow/format/notation across different types of algorithms

Different notations are required for different algorithms as the algorithms have been developed by many different people. This would require the painstaking process of reformatting the data to be compatible with each algorithm if multiple algorithms were tested.

2. Can easily modify preprocessing, algorithm choice, and hyper-parameter tuning making optimization easy

Modiying a piece of the overall process is now much easier than before because many of the steps are specified using the tidymodel packages in a convenient manner. Thus the entire process can be rerun after a simple change to preprocessing without much difficulty.

Splitting the Data

The first step after data exploration in machine learning analysis is to split the data into training and testing datasets.

The training dataset will be used to build and tune our model. This is the data that the model “learns” on.

The testing set will be used to evaluate the performance of our model in a more generalizable way. What do we mean by “generalizable”?

Remember that our main goal is to use our model to be able to predict air pollution levels in areas where there are no gravimetric monitors. Therefore, if our model is super good at predicting air pollution with the data that we use to build it, it might not do the best job for the areas where there are few to no monitors. This would cause us to have really good prediction accuracy and we might assume that we were going to do a good job estimating air pollution any time we use our model, but in fact this would likely not be the case. This situation is what we call overfitting .

Overfitting happens when we end up modeling not only the major relationships in our data but also the noise within our data.

knitr::include_graphics("https://miro.medium.com/max/1110/1*tBErXYVvTw2jSUYK7thU2A.png")

[source]

If we get fairly good prediction with our testing set then we will know that our model can be applied to other data and will perform fairly well. We will discuss this more later.

We will not touch the testing set until we have completed optimizing our model with the training set. This will allow us to have a less biased evaluation of how well our model can do with other data besides the data used in the training set to build the model. Ideally you would also want a completely independent dataset to further test the performance of your model.

Here is a great description of the differences between testing and training datasets.

We will use the rsample package to perform this step.

Theinitial_split() function allows us to specify how we want to split our data. Typically data is split into 3/4 for training and 1/4 for testing.This is the default proportion and does not need to be specified. However you can change the proportion using the prop argument, which we will do that here for illustrative purposes. You can also specify a variable to stratify by with the strata argument. This is useful if you have imbalanced categorical variables and you would like to intentionally make sure that there are similar number of samples of the rarer categories in both the testing and training sets. Otherwise the split is performed randomly.

The strata argument causes the random sampling to be conducted within the stratification variable. The can help ensure that the number of data points in the training data is equivalent to the proportions in the original data set.

In the case with our dataset, perhaps we would like our training set to have similar proportions of monitors from each of the states as in the initial data. This might be useful if we want our model to be generalizable across all of the states.

We can see that indeed there are different proportions of monitors in each state by using the count() function of the dpyr package.

# Scroll through the output!
count(pm, state) %>%
  print(n = 1e3)

# A tibble: 49 x 2
   state                    n
   <chr>                <int>
 1 Alabama                 24
 2 Arizona                 17
 3 Arkansas                16
 4 California              85
 5 Colorado                15
 6 Connecticut             14
 7 Delaware                 7
 8 District Of Columbia     3
 9 Florida                 29
10 Georgia                 28
11 Idaho                    7
12 Illinois                38
13 Indiana                 36
14 Iowa                    20
15 Kansas                  10
16 Kentucky                22
17 Louisiana               17
18 Maine                    1
19 Maryland                15
20 Massachusetts           16
21 Michigan                30
22 Minnesota               17
23 Mississippi             12
24 Missouri                13
25 Montana                 16
26 Nebraska                 7
27 Nevada                   4
28 New Hampshire            7
29 New Jersey              23
30 New Mexico              10
31 New York                24
32 North Carolina          35
33 North Dakota             4
34 Ohio                    44
35 Oklahoma                10
36 Oregon                  17
37 Pennsylvania            32
38 Rhode Island             5
39 South Carolina          14
40 South Dakota             9
41 Tennessee                3
42 Texas                   27
43 Utah                    14
44 Vermont                  4
45 Virginia                20
46 Washington               8
47 West Virginia           14
48 Wisconsin               21
49 Wyoming                 12

If our dataset were large enough it might be nice then to stratify by state, but our data is unfortunately not large enough. We will show how one would do this though for illustrative purposes. This option is often more important for classification applications of machine learning than it is for prediction.

Since the split is performed randomly, it is a good idea to use the set.seed() base function to ensure that if your rerun your code that your split will be the same next time. We can see the number of monitors in our training, testing, and original data by typing in the name of our split object. The result will look like this: <training data sample number, testing data sample number, original sample number>

set.seed(1234)
pm_split <-rsample::initial_split(data = pm, prop = 2/3)
pm_split

<Analysis/Assess/Total>
<584/292/876>

# If stratifying:
# pm_split_strata <-rsample::initial_split(data = pm, prop = 2/3, strata = "state")

Importantly the initial_split function only determines what rows of our pm data frame should be assigned for training or testing, it does not actually split the data.

To extract the testing and training data we can use the training() and testing() functions also of the rsample package.

 train_pm <-rsample::training(pm_split)
 test_pm <-rsample::testing(pm_split)
 
# Scroll through the output!
count(train_pm, state)

# A tibble: 48 x 2
   state                    n
   <chr>                <int>
 1 Alabama                 18
 2 Arizona                 12
 3 Arkansas                14
 4 California              54
 5 Colorado                12
 6 Connecticut              8
 7 Delaware                 6
 8 District Of Columbia     2
 9 Florida                 18
10 Georgia                 17
# … with 38 more rows

count(test_pm, state)

# A tibble: 48 x 2
   state                    n
   <chr>                <int>
 1 Alabama                  6
 2 Arizona                  5
 3 Arkansas                 2
 4 California              31
 5 Colorado                 3
 6 Connecticut              6
 7 Delaware                 1
 8 District Of Columbia     1
 9 Florida                 11
10 Georgia                 11
# … with 38 more rows

Variable Role Assignment and Preprocessing

In tidymodels we will create a recipe, which is a standardized format for a sequence of steps for processing the data.

This can be very useful because it makes testing out different pre-processing steps or different algorithms with the same pre-processing very easy and reproducible.

Creating a recipe specifies how a dataframe of predictors should be created - it specifies what variables to be used and the preprocessing steps but it does not execute these steps or create the dataframe of predictors.

List the ingredients / specify the variables with the recipe() function

The first thing to do to create a recipe is to specify which variables we will be using as our outcome and predictors using the recipe() function. In terms of the metaphor of baking, we can think of this as listing our ingredients. The naming convention for recipe object names is *_rec or rec.

In our case recall that our value variable, which is the average annual gravimetric monitor PM_2.5 concentration in ug/m³. Our predictors are all the other variables except the monitor ID, which is an id variable.

The reason not to include this variable is because this variable includes the county number and a number designating which particular monitor the values came from of the monitors there are in that county. Since this number is arbitrary and the county information is also given in the data, and the fact that each monitor only has one value in the value variable, nothing is gained by including this variable and it may instead introduce noise. However, it is useful to keep this data to take a look at what is happening later. We will show you what to do in this case in just a bit.

The simplest recipe with no preprocessing steps, would be to simply list the outcome and predictor variables.

We can do so in two ways:

1. Using formula notation
   
2. Assigning roles to each variable

Let’s look at the first way using formula notation, which looks like this:

outcome(s) ~ predictor(s)

If in the case of multiple predictors or a multivariate situation with two outcomes, use a plus sign

outcome1 + outcome2 ~ predictor1 + predictor2

If we want to include all predictors we can use a period like so:

outcome_variable_name ~ .

Now with our data we will start by making a recipe for our training data. In the simplest case we might use all predictors like this:

simple_rec <-train_pm %>%
  recipes::recipe(value ~ .)

simple_rec

Data Recipe

Inputs:

      role #variables
   outcome          1
 predictor         49

However, to deal with the id variable we could use the update_role() function of the recipes package. This option works well with the newer workflows package, however id variabels are often dropped from analyses that do not use this newer package as they can make the process difficult with using the parsnip package alone due to the fact that new levels (or possible values) may be introduced with the testing data.

simple_rec <-train_pm %>%
  recipes::recipe(value ~ .) %>%
  recipes::update_role(id, new_role = "id variable")

simple_rec

Data Recipe

Inputs:

        role #variables
 id variable          1
     outcome          1
   predictor         48

We could also specify the outcome and predictors in the same way as the id variable. Please see here for examples of other roles for variables. The role can be actually be any value.

The order is important here, as we first make all variables predictors and then override this role for the outcome and id variable. We will use the everything() function of the dplyr package to start with all of the variables in train_pm.

simple_rec <-recipe(train_pm) %>%
    update_role(everything(), new_role = "predictor")%>%
    update_role(value, new_role = "outcome")%>%
    update_role(id, new_role = "id variable")

simple_rec

Data Recipe

Inputs:

        role #variables
 id variable          1
     outcome          1
   predictor         48

If we want to take a look at our formula from our recipe we can do use the formula() function of the stats package.

formula(simple_rec)

value ~ fips + lat + lon + state + county + city + CMAQ + zcta + 
    zcta_area + zcta_pop + imp_a500 + imp_a1000 + imp_a5000 + 
    imp_a10000 + imp_a15000 + county_area + county_pop + log_dist_to_prisec + 
    log_pri_length_5000 + log_pri_length_10000 + log_pri_length_15000 + 
    log_pri_length_25000 + log_prisec_length_500 + log_prisec_length_1000 + 
    log_prisec_length_5000 + log_prisec_length_10000 + log_prisec_length_15000 + 
    log_prisec_length_25000 + log_nei_2008_pm25_sum_10000 + log_nei_2008_pm25_sum_15000 + 
    log_nei_2008_pm25_sum_25000 + log_nei_2008_pm10_sum_10000 + 
    log_nei_2008_pm10_sum_15000 + log_nei_2008_pm10_sum_25000 + 
    popdens_county + popdens_zcta + nohs + somehs + hs + somecollege + 
    associate + bachelor + grad + pov + hs_orless + urc2013 + 
    urc2006 + aod
<environment: 0x7fdfab277ab8>

We can also view our recipe in more detail using the base summary() function.

summary(simple_rec)

# A tibble: 50 x 4
   variable type    role        source  
   <chr>    <chr>   <chr>       <chr>   
 1 id       nominal id variable original
 2 value    numeric outcome     original
 3 fips     nominal predictor   original
 4 lat      numeric predictor   original
 5 lon      numeric predictor   original
 6 state    nominal predictor   original
 7 county   nominal predictor   original
 8 city     nominal predictor   original
 9 CMAQ     numeric predictor   original
10 zcta     nominal predictor   original
# … with 40 more rows

List the preprocessing steps using the step functions of the recipe package

The other thing the recipes package allows for is specifying pre-processing steps using a variety of step*() functions.

This link and this link show the many options for recipe step functions.

There are step functions for a variety of purposes:

1. Imputation – which means filling in missing values based on the existing data
2. Transformation – which means changing all values of a variable in the same way, typically to make it more normal or easier to interpret)
   
3. Discretization – which means converting continuous values into discrete or nominal values - binning for example to reduce the number of possible levels) (However this is generally not advisable!)
4. Encoding / Creating Dummy Variables – which means creating a numeric code for categorical variables More on Dummy Variables and one hot encoding
5. Data type conversions – which means changing from integer to factor or numeric to date etc.
6. Interaction term addition to the model – which means that we would be modeling for predictors that would influence the capacity of each other to predict the outcome
7. Normalization – which means centering and scaling the data to a similar range of values
8. Dimensionality Reduction/ Signal Extraction – which means mathematically obtaining a new smaller set of variables that capture the variation or signal in the original variables (ex. Principal Component Analysis and Independent Component Analysis)
9. Filtering – Filtering options for removing variables (ex. remove variables that are highly correlated to others or remove variables with very little variance and therefore likely little predictive capacity)
10. Row operations – which means performing functions on the values within the rows (ex. rearranging, filtering, imputing)
11. Checking functions – Sanity checks to look for missing values, to look at the variable classes etc.

All of the step functions look like step_* except for the check functions which look like check_*.

There are several ways to select what variables to apply steps to:
1) tidyselect methods: contains(), matches(), starts_with(), ends_with(), everything(), num_range()
2) based on the type: all_nominal(), all_numeric() , has_type() 3) based on the role: all_predictors(), all_outcomes(), has_role() 4) name - use the actual name of the variable/variables of interest

Let’s try adding some steps to our recipe.

We might consider log transforming our population and area variables (that aren’t densities) - let’s take a look at the range of these variables.

pm %>%
  select(matches("_pop|_area")) %>%
  map(range)

$zcta_area
[1]      15459 8164820625

$zcta_pop
[1]     0 95397

$county_area
[1]    33703512 51947229509

$county_pop
[1]     783 9818605

We can see that the range for each of these variables is quite large, we can log transform this data using the step_log() function of the recipes package.

We would also want to potentially one hot encode some of our categorical variables so that they can be used with certain algorithms. We can do this with the step_dummy() function and the one_hot = TRUE argument. Our fips variable includes a numeric code for state and county - and therefore is essentially a proxy for county. Since we already have county, we will just use it and keep the fips id as another ID variable.

We can remove the fips variable from the predictors using update_role() to make sure that the role is no longer "predictor". We can make the role anything we want actually, so we will keep it something identifiable.

We might also want to remove variables that appear to be redundant and are highly correlated with others, as we know from our exploratory data analysis that many of our variables are correlated with one another. We can do this using the step_corr() function.

It is also a good idea to remove variables with near-zero variance, which can be done with the step_nzv() function. Variables have low variance if all the values are very similar, the values are very sparse, or if they are highly imbalanced.

Examples where you might have near-zero variance variables include: ### I Struggled to find links for this… to back up my statements###

1. Similar Values - If the population density was nearly the same for every zcta that contained a monitor, then knowing the population density near our monitorwould contribute little to our model in assisting us to predict monitor air pollution values.
2. Sparse Data - If all of the monitors were in locations where the populations did not attend graduate school, then these values would mostly be zero, again this would do very little to help us distinguish our air pollution monitors.When many of the values are zero this is also called sparse data.
   
3. Imbalanced Data If nearly all of the monitors were located in one particular state, and all the others only had one monitor each, then the real predictive value would simply be in knowing if a monitor is located in that particular state or not. In this case we don’t want to remove our variable, we just want to simplify it.

It is important to add the steps to the recipe in an order that makes sense just like with a cooking recipe.

Thus first we are going to create numeric values for our categorical variables, then we will look at correlation and near-zero variance. We don’t want to remove some of our variables, like the CMAQ and aod variables so we can make sure they are kept in the model by excluding them from those steps. If we specifically wanted to remove a predictor we could use step_rm().

simple_rec %<>%
  update_role("fips", new_role = "county id") %>%
# create numeric dummy variables to encode for categorical variables
  step_dummy(state, county, city, zcta, one_hot = TRUE) %>%
  # identify and reomve all correlated predictors (now that they are numeric)
   step_corr(all_predictors(), - CMAQ, - aod)%>%
  # identify variables with near zero variance and remove
  step_nzv(all_predictors(), - CMAQ, - aod)
  
simple_rec

Data Recipe

Inputs:

        role #variables
   county id          1
 id variable          1
     outcome          1
   predictor         47

Operations:

Dummy variables from state, county, city, zcta
Correlation filter on all_predictors, -, CMAQ, -, aod
Sparse, unbalanced variable filter on all_predictors, -, CMAQ, -, aod

Running the preprocessing

The next major function of the recipes package is prep().

This function updates the recipe object based on the training data. It estimates parameters (estimating the required quantities and statistics required by the steps for the variables) for preprocessing and updates the model terms, as some of the predictors may be removed, this allows the recipe to be ready to use on other datasets. It doesn’t necessarily actually execute the preprocessing itself, however we will specify in argument for it to do this so that we can take a look at the preprocessed data.

There are some important arguments to know about: 1) training - you must supply a training data set to estimate parameters for preprocessing operations (recipe steps) - this may already be included in your recipe - as is the case for us 2) fresh - if TRUE - will retrain and estimate parameters for any previous steps that were already prepped if you add more steps to the recipe 3) verbose - if TRUE shows the progress as the steps are evaluated and the size of the preprocessed training set 4) retain - if TRUE then the preprocessed training set will be saved within the recipe (as template). This is good if you are likely to add more steps and don’t want to rerun the prep() on the previous steps. However this can make the recipe size large. This is necessary if you want to actually look at the preprocessed data.

prepped_rec <- prep(simple_rec, verbose = TRUE, retain = TRUE )

oper 1 step dummy [training] 
oper 2 step corr [training] 
oper 3 step nzv [training] 
The retained training set is ~ 0.26 Mb  in memory.

names(prepped_rec)

[1] "var_info"       "term_info"      "steps"          "template"      
[5] "levels"         "retained"       "tr_info"        "orig_lvls"     
[9] "last_term_info"

There are also lots of useful things to checkout in the output of prep(). You can see: 1) the steps that were run
2) the variable info (var_info)
3) the model term_info 4) the new levels of the variables 5) the original levels of the variables orig_lvls
6) info about the training data set size and completeness (tr_info)

Note: You may see the prep.recipe() function in material that you read about the recipes package. This is referring to the prep() function of the recipes package.

Extracting the preprocessed training data

Since we retained our preprocessed training data, we can take a look at it like by using the juice() function of the recipes package like this:

# Scroll therough the output!
juiced_train<- juice(prepped_rec)
glimpse(juiced_train)

Rows: 584
Columns: 36
$ id                          <fct> 1003.001, 1027.0001, 1033.1002, 1055.001,…
$ value                       <dbl> 9.597647, 10.800000, 11.212174, 12.375394…
$ fips                        <fct> 1003, 1027, 1033, 1055, 1069, 1073, 1073,…
$ lat                         <dbl> 30.49800, 33.28126, 34.75878, 33.99375, 3…
$ lon                         <dbl> -87.88141, -85.80218, -87.65056, -85.9910…
$ CMAQ                        <dbl> 8.098836, 9.766208, 9.402679, 9.241744, 9…
$ zcta_area                   <dbl> 190980522, 374132430, 16716984, 154069359…
$ zcta_pop                    <dbl> 27829, 5103, 9042, 20045, 30217, 9010, 16…
$ imp_a500                    <dbl> 0.01730104, 1.96972318, 19.17301038, 16.4…
$ imp_a15000                  <dbl> 1.4386207, 0.3359198, 5.2472094, 5.161210…
$ county_area                 <dbl> 4117521611, 1564252280, 1534877333, 13856…
$ county_pop                  <dbl> 182265, 13932, 54428, 104430, 101547, 658…
$ log_dist_to_prisec          <dbl> 4.648181, 7.219907, 5.760131, 5.261457, 7…
$ log_pri_length_5000         <dbl> 8.517193, 8.517193, 8.517193, 9.066563, 8…
$ log_pri_length_25000        <dbl> 11.32735, 10.12663, 10.15769, 12.01356, 1…
$ log_prisec_length_500       <dbl> 7.295356, 6.214608, 8.611945, 8.740680, 6…
$ log_prisec_length_1000      <dbl> 8.195119, 7.600902, 9.735569, 9.627898, 7…
$ log_prisec_length_5000      <dbl> 10.815042, 10.170878, 11.770407, 11.72888…
$ log_prisec_length_10000     <dbl> 11.886803, 11.405543, 12.840663, 12.76827…
$ log_nei_2008_pm10_sum_15000 <dbl> 2.26783411, 3.31111648, 6.70127741, 4.462…
$ log_nei_2008_pm10_sum_25000 <dbl> 5.628728, 3.311116, 7.148858, 4.678311, 3…
$ popdens_county              <dbl> 44.265706, 8.906492, 35.460814, 75.367038…
$ popdens_zcta                <dbl> 145.7164307, 13.6395554, 540.8870404, 130…
$ nohs                        <dbl> 3.3, 11.6, 7.3, 4.3, 5.8, 7.1, 2.7, 11.1,…
$ somehs                      <dbl> 4.9, 19.1, 15.8, 13.3, 11.6, 17.1, 6.6, 1…
$ hs                          <dbl> 25.1, 33.9, 30.6, 27.8, 29.8, 37.2, 30.7,…
$ somecollege                 <dbl> 19.7, 18.8, 20.9, 29.2, 21.4, 23.5, 25.7,…
$ associate                   <dbl> 8.2, 8.0, 7.6, 10.1, 7.9, 7.3, 8.0, 4.1, …
$ bachelor                    <dbl> 25.3, 5.5, 12.7, 10.0, 13.7, 5.9, 17.6, 7…
$ grad                        <dbl> 13.5, 3.1, 5.1, 5.4, 9.8, 2.0, 8.7, 2.9, …
$ pov                         <dbl> 6.1, 19.5, 19.0, 8.8, 15.6, 25.5, 7.3, 8.…
$ hs_orless                   <dbl> 33.3, 64.6, 53.7, 45.4, 47.2, 61.4, 40.0,…
$ urc2013                     <dbl> 4, 6, 4, 4, 4, 1, 1, 1, 1, 1, 2, 3, 3, 3,…
$ aod                         <dbl> 37.363636, 34.818182, 36.000000, 43.41666…
$ state_California            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ city_Not.in.a.city          <dbl> 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0,…

For easy comparison sake - here is our original data: #### {.scrollable }

# Scroll therough the output!
glimpse(pm)

Rows: 876
Columns: 50
$ id                          <fct> 1003.001, 1027.0001, 1033.1002, 1049.1003…
$ value                       <dbl> 9.597647, 10.800000, 11.212174, 11.659091…
$ fips                        <fct> 1003, 1027, 1033, 1049, 1055, 1069, 1073,…
$ lat                         <dbl> 30.49800, 33.28126, 34.75878, 34.28763, 3…
$ lon                         <dbl> -87.88141, -85.80218, -87.65056, -85.9683…
$ state                       <chr> "Alabama", "Alabama", "Alabama", "Alabama…
$ county                      <chr> "Baldwin", "Clay", "Colbert", "DeKalb", "…
$ city                        <chr> "Fairhope", "Ashland", "Muscle Shoals", "…
$ CMAQ                        <dbl> 8.098836, 9.766208, 9.402679, 8.534772, 9…
$ zcta                        <fct> 36532, 36251, 35660, 35962, 35901, 36303,…
$ zcta_area                   <dbl> 190980522, 374132430, 16716984, 203836235…
$ zcta_pop                    <dbl> 27829, 5103, 9042, 8300, 20045, 30217, 90…
$ imp_a500                    <dbl> 0.01730104, 1.96972318, 19.17301038, 5.78…
$ imp_a1000                   <dbl> 1.4096021, 0.8531574, 11.1448962, 3.86764…
$ imp_a5000                   <dbl> 3.3360118, 0.9851479, 15.1786154, 1.23114…
$ imp_a10000                  <dbl> 1.9879187, 0.5208189, 9.7253870, 1.031646…
$ imp_a15000                  <dbl> 1.4386207, 0.3359198, 5.2472094, 0.973044…
$ county_area                 <dbl> 4117521611, 1564252280, 1534877333, 20126…
$ county_pop                  <dbl> 182265, 13932, 54428, 71109, 104430, 1015…
$ log_dist_to_prisec          <dbl> 4.648181, 7.219907, 5.760131, 3.721489, 5…
$ log_pri_length_5000         <dbl> 8.517193, 8.517193, 8.517193, 8.517193, 9…
$ log_pri_length_10000        <dbl> 9.210340, 9.210340, 9.274303, 10.409411, …
$ log_pri_length_15000        <dbl> 9.630228, 9.615805, 9.658899, 11.173626, …
$ log_pri_length_25000        <dbl> 11.32735, 10.12663, 10.15769, 11.90959, 1…
$ log_prisec_length_500       <dbl> 7.295356, 6.214608, 8.611945, 7.310155, 8…
$ log_prisec_length_1000      <dbl> 8.195119, 7.600902, 9.735569, 8.585843, 9…
$ log_prisec_length_5000      <dbl> 10.815042, 10.170878, 11.770407, 10.21420…
$ log_prisec_length_10000     <dbl> 11.88680, 11.40554, 12.84066, 11.50894, 1…
$ log_prisec_length_15000     <dbl> 12.205723, 12.042963, 13.282656, 12.35366…
$ log_prisec_length_25000     <dbl> 13.41395, 12.79980, 13.79973, 13.55979, 1…
$ log_nei_2008_pm25_sum_10000 <dbl> 0.318035438, 3.218632928, 6.573127301, 0.…
$ log_nei_2008_pm25_sum_15000 <dbl> 1.967358961, 3.218632928, 6.581917457, 3.…
$ log_nei_2008_pm25_sum_25000 <dbl> 5.067308, 3.218633, 6.875900, 4.887665, 4…
$ log_nei_2008_pm10_sum_10000 <dbl> 1.35588511, 3.31111648, 6.69187313, 0.000…
$ log_nei_2008_pm10_sum_15000 <dbl> 2.26783411, 3.31111648, 6.70127741, 3.350…
$ log_nei_2008_pm10_sum_25000 <dbl> 5.628728, 3.311116, 7.148858, 5.171920, 4…
$ popdens_county              <dbl> 44.265706, 8.906492, 35.460814, 35.330814…
$ popdens_zcta                <dbl> 145.716431, 13.639555, 540.887040, 40.718…
$ nohs                        <dbl> 3.3, 11.6, 7.3, 14.3, 4.3, 5.8, 7.1, 2.7,…
$ somehs                      <dbl> 4.9, 19.1, 15.8, 16.7, 13.3, 11.6, 17.1, …
$ hs                          <dbl> 25.1, 33.9, 30.6, 35.0, 27.8, 29.8, 37.2,…
$ somecollege                 <dbl> 19.7, 18.8, 20.9, 14.9, 29.2, 21.4, 23.5,…
$ associate                   <dbl> 8.2, 8.0, 7.6, 5.5, 10.1, 7.9, 7.3, 8.0, …
$ bachelor                    <dbl> 25.3, 5.5, 12.7, 7.9, 10.0, 13.7, 5.9, 17…
$ grad                        <dbl> 13.5, 3.1, 5.1, 5.8, 5.4, 9.8, 2.0, 8.7, …
$ pov                         <dbl> 6.1, 19.5, 19.0, 13.8, 8.8, 15.6, 25.5, 7…
$ hs_orless                   <dbl> 33.3, 64.6, 53.7, 66.0, 45.4, 47.2, 61.4,…
$ urc2013                     <dbl> 4, 6, 4, 6, 4, 4, 1, 1, 1, 1, 1, 1, 1, 2,…
$ urc2006                     <dbl> 5, 6, 4, 5, 4, 4, 1, 1, 1, 1, 1, 1, 1, 2,…
$ aod                         <dbl> 37.36364, 34.81818, 36.00000, 33.08333, 4…

Notice how we only have 36 variables now instead of 50! Two of these are our ID variables (fips and the actual monitor ID (id)) and one is our outcome (value). Thus we only have 33 predictors now. We can also see that variables that we no longer have any categorical variables. Variables like state are gone and only state_California remains as it was the only state identity to have nonzero variance. We can see that California had the largest number of monitors compared to the other states. We can also see that there were more monitors listed as "Not in a city" than any city.

pm %>% count(state)

# A tibble: 49 x 2
   state                    n
   <chr>                <int>
 1 Alabama                 24
 2 Arizona                 17
 3 Arkansas                16
 4 California              85
 5 Colorado                15
 6 Connecticut             14
 7 Delaware                 7
 8 District Of Columbia     3
 9 Florida                 29
10 Georgia                 28
# … with 39 more rows

pm %>% count(city)

# A tibble: 607 x 2
   city            n
   <chr>       <int>
 1 Aberdeen        1
 2 Akron           2
 3 Albany          3
 4 Albuquerque     2
 5 Alexandria      1
 6 Allen Park      1
 7 Altamont        1
 8 Alton           1
 9 Amarillo        1
10 Anadarko        1
# … with 597 more rows

Note: Recall that you must specify retain = TRUE argument of the prep() function to use juice().

Extracting the preprocessed testing data

Accroding to the tidymodels documentation:

bake() takes a trained recipe and applies the operations to a data set to create a design matrix. for example: it applies the centering to new data sets using these means used to create the recipe

If you wanted to look at the preprocessed testing data you would use the bake() function of the recipes package. (You generally want to leave your testing data alone, but it is good to look for issues like the introduction of NA values).

# Scroll therough the output!
baked_test_pm <- recipes::bake(prepped_rec, new_data = test_pm)
glimpse(baked_test_pm )

Rows: 292
Columns: 36
$ id                          <fct> 1049.1003, 1073.101, 1073.2006, 1089.0014…
$ value                       <dbl> 11.659091, 13.114545, 12.228125, 12.23294…
$ fips                        <fct> 1049, 1073, 1073, 1089, 1103, 1121, 4013,…
$ lat                         <dbl> 34.28763, 33.54528, 33.38639, 34.68767, 3…
$ lon                         <dbl> -85.96830, -86.54917, -86.81667, -86.5863…
$ CMAQ                        <dbl> 8.534772, 9.303766, 10.235612, 9.343611, …
$ zcta_area                   <dbl> 203836235, 148994881, 56063756, 46963946,…
$ zcta_pop                    <dbl> 8300, 14212, 32390, 21297, 30545, 7713, 5…
$ imp_a500                    <dbl> 5.78200692, 0.06055363, 42.42820069, 23.2…
$ imp_a15000                  <dbl> 0.9730444, 2.9956557, 12.7487614, 10.3555…
$ county_area                 <dbl> 2012662359, 2878192209, 2878192209, 20761…
$ county_pop                  <dbl> 71109, 658466, 658466, 334811, 119490, 82…
$ log_dist_to_prisec          <dbl> 3.721489, 7.301545, 4.721755, 4.659519, 6…
$ log_pri_length_5000         <dbl> 8.517193, 9.683336, 10.737240, 8.517193, …
$ log_pri_length_25000        <dbl> 11.90959, 12.53777, 12.99669, 11.47391, 1…
$ log_prisec_length_500       <dbl> 7.310155, 6.214608, 7.528913, 8.760549, 6…
$ log_prisec_length_1000      <dbl> 8.585843, 7.600902, 9.342290, 9.543183, 8…
$ log_prisec_length_5000      <dbl> 10.214200, 11.262645, 11.713190, 11.48606…
$ log_prisec_length_10000     <dbl> 11.50894, 12.14101, 12.53899, 12.68440, 1…
$ log_nei_2008_pm10_sum_15000 <dbl> 3.3500444, 6.6241114, 5.8268686, 3.861625…
$ log_nei_2008_pm10_sum_25000 <dbl> 5.1719202, 7.5490587, 8.8205542, 5.219092…
$ popdens_county              <dbl> 35.330814, 228.777633, 228.777633, 161.26…
$ popdens_zcta                <dbl> 40.718962, 95.385827, 577.735106, 453.475…
$ nohs                        <dbl> 14.3, 7.2, 0.8, 1.2, 4.8, 16.7, 19.1, 6.4…
$ somehs                      <dbl> 16.7, 12.2, 2.6, 3.1, 7.8, 33.3, 15.6, 9.…
$ hs                          <dbl> 35.0, 32.2, 12.9, 15.1, 28.7, 37.5, 26.5,…
$ somecollege                 <dbl> 14.9, 19.0, 17.9, 20.5, 25.0, 12.5, 18.0,…
$ associate                   <dbl> 5.5, 6.8, 5.2, 6.5, 7.5, 0.0, 6.0, 8.8, 3…
$ bachelor                    <dbl> 7.9, 14.8, 35.5, 30.4, 18.2, 0.0, 10.6, 1…
$ grad                        <dbl> 5.8, 7.7, 25.2, 23.3, 8.0, 0.0, 4.1, 5.7,…
$ pov                         <dbl> 13.8, 10.5, 2.1, 5.2, 8.3, 18.8, 21.4, 14…
$ hs_orless                   <dbl> 66.0, 51.6, 16.3, 19.4, 41.3, 87.5, 61.2,…
$ urc2013                     <dbl> 6, 1, 1, 3, 4, 5, 1, 2, 5, 4, 4, 6, 6, 1,…
$ aod                         <dbl> 33.08333, 42.45455, 44.25000, 42.41667, 4…
$ state_California            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,…
$ city_Not.in.a.city          <dbl> NA, NA, NA, NA, NA, NA, 0, NA, NA, NA, NA…

Notice that our city_Not.in.a.city variable seems to be NA values. Why might that be?

Ah! Perphas it is because some of our levels were not previously seen in the training set!

Let’s take a look using the set operations of the dplyr package. We can take a look at cities that were different between the test and training set.

traincities <- train_pm %>% distinct(city)
testcities <- test_pm %>% distinct(city)

#get the number of cities that were different
dim(dplyr::setdiff(traincities, testcities))

[1] 376   1

#get the number of cities that overlapped
dim(dplyr::intersect(traincities, testcities))

[1] 51  1

Indeed, there are lots of different cities in our test data that are not in our training data!

Maybe remove this?: Thus we need to update our original recipe to include a very important step function called step_novel() this helps in cases like this were there are new factors in our testing set that were not in our training set. It is a good idea to include this in most of your recipes where you have a categorical variables with many distinct values. This step needs to come before we create dummy variables. However, we are also creating a dummy variable from this, which still results in a problem.

Let’s modify the city variable to be values of in a city or not in a city using the if_else() function of dplyr. Alternatively you could create a custom step function to do this and add the step function to your recipe, but that is beyond the scope of this case study.

We need to create a new recipe to move forward, as the levels of our variables are established then. We would also potentially have this issue for state and county. So let’s also do a similar thing for state. The county variables appears to get dropped due to either correlation or near zero variance. It is likely due to near zero variance becuase this is the more granular of these geographic categorical variables and likely sparse.

pm[["city"]]<-if_else(
  pm[["city"]] == "Not in a city", 
  true = "Not in a city", false = "In a city")

pm[["state"]] <- if_else(
  pm[["state"]] == "California", 
  true = "California", false = "Not California")

set.seed(1234) # same seed as before

pm_split <-rsample::initial_split(data = pm, prop = 2/3)
pm_split

<Analysis/Assess/Total>
<584/292/876>

 train_pm <-rsample::training(pm_split)
 test_pm <-rsample::testing(pm_split)

  
novel_rec <-recipe(train_pm) %>%
    update_role(everything(), new_role = "predictor")%>%
    update_role(value, new_role = "outcome")%>%
    update_role(id, new_role = "id variable")%>%
    update_role("fips", new_role = "county id")%>%
# create numeric dummy variables to encode for categorical variables
    step_dummy(state, county, city, zcta, one_hot = TRUE) %>%
# identify and reomve all correlated predictors (now that they are numeric)
   step_corr(all_numeric())%>%
# identify variables with near zero variance and remove
  step_nzv(all_numeric()) # 

Now let’s retrain our training data and try baking our test data:

prepped_rec <- prep(novel_rec, verbose = TRUE, retain = TRUE)

oper 1 step dummy [training] 
oper 2 step corr [training] 
oper 3 step nzv [training] 
The retained training set is ~ 0.26 Mb  in memory.

# Scroll therough the output!
juiced_train<- juice(prepped_rec)
glimpse(juiced_train)

Rows: 584
Columns: 37
$ id                          <fct> 1003.001, 1027.0001, 1033.1002, 1055.001,…
$ value                       <dbl> 9.597647, 10.800000, 11.212174, 12.375394…
$ fips                        <fct> 1003, 1027, 1033, 1055, 1069, 1073, 1073,…
$ lat                         <dbl> 30.49800, 33.28126, 34.75878, 33.99375, 3…
$ lon                         <dbl> -87.88141, -85.80218, -87.65056, -85.9910…
$ CMAQ                        <dbl> 8.098836, 9.766208, 9.402679, 9.241744, 9…
$ zcta_area                   <dbl> 190980522, 374132430, 16716984, 154069359…
$ zcta_pop                    <dbl> 27829, 5103, 9042, 20045, 30217, 9010, 16…
$ imp_a500                    <dbl> 0.01730104, 1.96972318, 19.17301038, 16.4…
$ imp_a15000                  <dbl> 1.4386207, 0.3359198, 5.2472094, 5.161210…
$ county_area                 <dbl> 4117521611, 1564252280, 1534877333, 13856…
$ county_pop                  <dbl> 182265, 13932, 54428, 104430, 101547, 658…
$ log_dist_to_prisec          <dbl> 4.648181, 7.219907, 5.760131, 5.261457, 7…
$ log_pri_length_5000         <dbl> 8.517193, 8.517193, 8.517193, 9.066563, 8…
$ log_pri_length_25000        <dbl> 11.32735, 10.12663, 10.15769, 12.01356, 1…
$ log_prisec_length_500       <dbl> 7.295356, 6.214608, 8.611945, 8.740680, 6…
$ log_prisec_length_1000      <dbl> 8.195119, 7.600902, 9.735569, 9.627898, 7…
$ log_prisec_length_5000      <dbl> 10.815042, 10.170878, 11.770407, 11.72888…
$ log_prisec_length_10000     <dbl> 11.886803, 11.405543, 12.840663, 12.76827…
$ log_prisec_length_25000     <dbl> 13.41395, 12.79980, 13.79973, 13.70026, 1…
$ log_nei_2008_pm10_sum_15000 <dbl> 2.26783411, 3.31111648, 6.70127741, 4.462…
$ log_nei_2008_pm10_sum_25000 <dbl> 5.628728, 3.311116, 7.148858, 4.678311, 3…
$ popdens_county              <dbl> 44.265706, 8.906492, 35.460814, 75.367038…
$ popdens_zcta                <dbl> 145.7164307, 13.6395554, 540.8870404, 130…
$ nohs                        <dbl> 3.3, 11.6, 7.3, 4.3, 5.8, 7.1, 2.7, 11.1,…
$ somehs                      <dbl> 4.9, 19.1, 15.8, 13.3, 11.6, 17.1, 6.6, 1…
$ hs                          <dbl> 25.1, 33.9, 30.6, 27.8, 29.8, 37.2, 30.7,…
$ somecollege                 <dbl> 19.7, 18.8, 20.9, 29.2, 21.4, 23.5, 25.7,…
$ associate                   <dbl> 8.2, 8.0, 7.6, 10.1, 7.9, 7.3, 8.0, 4.1, …
$ bachelor                    <dbl> 25.3, 5.5, 12.7, 10.0, 13.7, 5.9, 17.6, 7…
$ grad                        <dbl> 13.5, 3.1, 5.1, 5.4, 9.8, 2.0, 8.7, 2.9, …
$ pov                         <dbl> 6.1, 19.5, 19.0, 8.8, 15.6, 25.5, 7.3, 8.…
$ hs_orless                   <dbl> 33.3, 64.6, 53.7, 45.4, 47.2, 61.4, 40.0,…
$ urc2013                     <dbl> 4, 6, 4, 4, 4, 1, 1, 1, 1, 1, 2, 3, 3, 3,…
$ aod                         <dbl> 37.363636, 34.818182, 36.000000, 43.41666…
$ state_Not.California        <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ city_Not.in.a.city          <dbl> 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0,…

Notice, it looks like we gained the log_prisec_length_25000 back with this recipe using the data with our changes to state and city.

baked_test_pm<- recipes::bake(prepped_rec, new_data = test_pm)
glimpse(baked_test_pm)

Rows: 292
Columns: 37
$ id                          <fct> 1049.1003, 1073.101, 1073.2006, 1089.0014…
$ value                       <dbl> 11.659091, 13.114545, 12.228125, 12.23294…
$ fips                        <fct> 1049, 1073, 1073, 1089, 1103, 1121, 4013,…
$ lat                         <dbl> 34.28763, 33.54528, 33.38639, 34.68767, 3…
$ lon                         <dbl> -85.96830, -86.54917, -86.81667, -86.5863…
$ CMAQ                        <dbl> 8.534772, 9.303766, 10.235612, 9.343611, …
$ zcta_area                   <dbl> 203836235, 148994881, 56063756, 46963946,…
$ zcta_pop                    <dbl> 8300, 14212, 32390, 21297, 30545, 7713, 5…
$ imp_a500                    <dbl> 5.78200692, 0.06055363, 42.42820069, 23.2…
$ imp_a15000                  <dbl> 0.9730444, 2.9956557, 12.7487614, 10.3555…
$ county_area                 <dbl> 2012662359, 2878192209, 2878192209, 20761…
$ county_pop                  <dbl> 71109, 658466, 658466, 334811, 119490, 82…
$ log_dist_to_prisec          <dbl> 3.721489, 7.301545, 4.721755, 4.659519, 6…
$ log_pri_length_5000         <dbl> 8.517193, 9.683336, 10.737240, 8.517193, …
$ log_pri_length_25000        <dbl> 11.90959, 12.53777, 12.99669, 11.47391, 1…
$ log_prisec_length_500       <dbl> 7.310155, 6.214608, 7.528913, 8.760549, 6…
$ log_prisec_length_1000      <dbl> 8.585843, 7.600902, 9.342290, 9.543183, 8…
$ log_prisec_length_5000      <dbl> 10.214200, 11.262645, 11.713190, 11.48606…
$ log_prisec_length_10000     <dbl> 11.50894, 12.14101, 12.53899, 12.68440, 1…
$ log_prisec_length_25000     <dbl> 13.55979, 14.08915, 14.27363, 13.87170, 1…
$ log_nei_2008_pm10_sum_15000 <dbl> 3.3500444, 6.6241114, 5.8268686, 3.861625…
$ log_nei_2008_pm10_sum_25000 <dbl> 5.1719202, 7.5490587, 8.8205542, 5.219092…
$ popdens_county              <dbl> 35.330814, 228.777633, 228.777633, 161.26…
$ popdens_zcta                <dbl> 40.718962, 95.385827, 577.735106, 453.475…
$ nohs                        <dbl> 14.3, 7.2, 0.8, 1.2, 4.8, 16.7, 19.1, 6.4…
$ somehs                      <dbl> 16.7, 12.2, 2.6, 3.1, 7.8, 33.3, 15.6, 9.…
$ hs                          <dbl> 35.0, 32.2, 12.9, 15.1, 28.7, 37.5, 26.5,…
$ somecollege                 <dbl> 14.9, 19.0, 17.9, 20.5, 25.0, 12.5, 18.0,…
$ associate                   <dbl> 5.5, 6.8, 5.2, 6.5, 7.5, 0.0, 6.0, 8.8, 3…
$ bachelor                    <dbl> 7.9, 14.8, 35.5, 30.4, 18.2, 0.0, 10.6, 1…
$ grad                        <dbl> 5.8, 7.7, 25.2, 23.3, 8.0, 0.0, 4.1, 5.7,…
$ pov                         <dbl> 13.8, 10.5, 2.1, 5.2, 8.3, 18.8, 21.4, 14…
$ hs_orless                   <dbl> 66.0, 51.6, 16.3, 19.4, 41.3, 87.5, 61.2,…
$ urc2013                     <dbl> 6, 1, 1, 3, 4, 5, 1, 2, 5, 4, 4, 6, 6, 1,…
$ aod                         <dbl> 33.08333, 42.45455, 44.25000, 42.41667, 4…
$ state_Not.California        <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,…
$ city_Not.in.a.city          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…

Great now we no longer have NA values! :)

Note: if you use the skip option for some of the preprocessing steps, be careful. juice() will show all of the results ignoring skip = TRUE. bake() will not necessarily conduct these steps on the new data.

Specifying the Model

So far we have used rsample to split the data and recipes to assign variable and to specify and prep our preprocessing (as well as to optionally extract the preprocessed data).

We will now use the parsnip package (which is similar to the previous caret package - and hence why it is named after the vegetable) to specify our model.

There are four aspects to define about our model:
1) the type of model (using specific funtions in parsnip like rand_forest(), logistic_reg() etc.)
2) the mode of learning - classification or regression (using the set_mode() function)
3) the package or engine that we will use to implement the type of model selected (using the set_engine() function)
4) any arguments necessary for the model/package selected (using the set_args()function - for example the mtry = argument for random forest which is the number of variables to be used as options for splitting at each tree node)

We are going to start our analysis with a linear regression but we will demonstrate how we can try different models.

The first thing we do is define what type of model we would like to use. See here for modeling [options]in parsnip.

PM_model <- parsnip::linear_reg() #PM for particulate mater
PM_model

Linear Regression Model Specification (regression)

OK. So far all we have told parsnip is we want to use a linear regression… Let’s tell parsnip more about what we want.

We would like to use the ordinary least squares method to fit our linear regression. So we will tell parsnip that we want to use the lm package to implement our linear regression (there are many options actually- such as rstan glmnet, keras, and sparklyr). We will do so by using the set_engine() function of the parsnip package.

lm_PM_model <- 
  PM_model  %>%
  parsnip::set_engine("lm")

lm_PM_model

Linear Regression Model Specification (regression)

Computational engine: lm 

In some cases some packages can do either classification or prediction, so it is a good idea to specify which mode you intend to perform. You can do this with the set_mode() function of the parsnip package, by using either set_mode("classification") or set_mode("regression").

lm_PM_model <- 
  PM_model  %>%
  parsnip::set_engine("lm") %>%
  set_mode("regression")

lm_PM_model

Linear Regression Model Specification (regression)

Computational engine: lm 

Fitting the Model: two ways - workflows and parsnip

To fit our model we can use the parsnip package and then assess our fit using the yardstick package.

However a newer package called workflows allows us to keep track of both our preprocessing steps and our model specification. It also allows us to implement fancier optimizations in an automated way and it is currently being developed to also handle post-processing operations, so it is good to learn about it!

So we will now create a workflow with the recipe (our preprocessing specifications) that we made and the model that we just specified.

First we use the workflow() function of the workflows package to create a workflow.

Then we add our recipe with the add_recipe() function and we add our model with the add_model() function of the workflows package.

Note: We do not need to actually prep our recipe before using workflows!

PM_wflow <-workflows::workflow() %>%
           workflows::add_recipe(novel_rec) %>%
           workflows::add_model(lm_PM_model)
PM_wflow

══ Workflow ═══════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: linear_reg()

── Preprocessor ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
3 Recipe Steps

● step_dummy()
● step_corr()
● step_nzv()

── Model ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Linear Regression Model Specification (regression)

Computational engine: lm 

Ah, nice. Notice how it tells us about both our preprocessing steps and our model specifications.

Now we can prepare the recipe (estimate the parameters) and fit the model to our training data all at once. Printing the output we can see the coefficients of the model.

PM_wflow_fit <- parsnip::fit(PM_wflow, data = train_pm)
PM_wflow_fit

══ Workflow [trained] ═════════════════════════════════════════════════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: linear_reg()

── Preprocessor ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
3 Recipe Steps

● step_dummy()
● step_corr()
● step_nzv()

── Model ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

Call:
stats::lm(formula = formula, data = data)

Coefficients:
                (Intercept)                          lat  
                  2.265e+02                    3.639e-02  
                        lon                         CMAQ  
                  2.579e-02                    2.847e-01  
                  zcta_area                     zcta_pop  
                  3.638e-10                    7.880e-06  
                   imp_a500                   imp_a15000  
                  7.453e-03                   -1.140e-03  
                county_area                   county_pop  
                 -2.116e-11                   -2.156e-07  
         log_dist_to_prisec          log_pri_length_5000  
                  4.064e-02                   -1.842e-01  
       log_pri_length_25000        log_prisec_length_500  
                  7.936e-03                    2.572e-01  
     log_prisec_length_1000       log_prisec_length_5000  
                 -2.288e-02                    4.743e-01  
    log_prisec_length_10000      log_prisec_length_25000  
                 -1.410e-01                    4.658e-01  
log_nei_2008_pm10_sum_15000  log_nei_2008_pm10_sum_25000  
                  1.088e-01                    5.255e-02  
             popdens_county                 popdens_zcta  
                 -4.114e-05                   -1.816e-05  
                       nohs                       somehs  
                 -2.234e+00                   -2.270e+00  
                         hs                  somecollege  
                 -2.271e+00                   -2.274e+00  
                  associate                     bachelor  
                 -2.272e+00                   -2.277e+00  
                       grad                          pov  
                 -2.284e+00                    7.001e-03  
                  hs_orless                      urc2013  
                         NA                    2.153e-01  
                        aod         state_Not.California  
                  2.555e-02                   -3.445e+00  
         city_Not.in.a.city  
                  3.550e-01  

Otherwise we could have done this without the workflows package. Notice here we will used the processed training data (juiced_train) as opposed to the raw training data that we used with the workflow we created with workflows.

In this case, we actually need to write your model again! Recall that id and fips are ID variables and that values is our outcome of interest (the pm air pollution measure at each monitor).

juiced_train_ready <- juiced_train%>% select( -id, -fips)
PM_fit <- lm_PM_model %>% 
parsnip::fit(value ~., data =juiced_train_ready)

Looking at model fit with broom

The broom package allows for an easy/tidy way to look at the fitted model:

tidy() grabs the coefficients from the model
glance() summarizes the model fit and gives us an idea about how well the model might perform augment() gives a 150 row observation level summary of the data and fit

These broom functions currently only work with parsnip objects not raw workflows objects. To use the tidy() function with workflows we need to first use the pull_workflow_fit() function.

broom::tidy(PM_fit) %>% arrange(p.value)

# A tibble: 34 x 5
   term                     estimate    std.error statistic  p.value
   <chr>                       <dbl>        <dbl>     <dbl>    <dbl>
 1 state_Not.California -3.44        0.436            -7.91 1.45e-14
 2 CMAQ                  0.285       0.0430            6.61 8.90e-11
 3 aod                   0.0256      0.00575           4.45 1.06e- 5
 4 lon                   0.0258      0.00998           2.58 1.00e- 2
 5 county_pop           -0.000000216 0.0000000934     -2.31 2.14e- 2
 6 urc2013               0.215       0.101             2.13 3.35e- 2
 7 grad                 -2.28        1.20             -1.90 5.78e- 2
 8 bachelor             -2.28        1.20             -1.90 5.84e- 2
 9 somecollege          -2.27        1.20             -1.89 5.87e- 2
10 associate            -2.27        1.20             -1.89 5.90e- 2
# … with 24 more rows

broom::glance(PM_fit[["fit"]])

# A tibble: 1 x 11
  r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>    <dbl> <int>  <dbl> <dbl> <dbl>
1     0.433         0.399  2.08      12.7 1.02e-48    34 -1240. 2550. 2703.
# … with 2 more variables: deviance <dbl>, df.residual <int>

broom::augment(PM_fit[["fit"]]) # this also gives us the fitted values, standard error for each and more!

# A tibble: 584 x 42
   value   lat   lon  CMAQ zcta_area zcta_pop imp_a500 imp_a15000 county_area
   <dbl> <dbl> <dbl> <dbl>     <dbl>    <dbl>    <dbl>      <dbl>       <dbl>
 1  9.60  30.5 -87.9  8.10 190980522    27829   0.0173      1.44   4117521611
 2 10.8   33.3 -85.8  9.77 374132430     5103   1.97        0.336  1564252280
 3 11.2   34.8 -87.7  9.40  16716984     9042  19.2         5.25   1534877333
 4 12.4   34.0 -86.0  9.24 154069359    20045  16.5         5.16   1385618994
 5 10.5   31.2 -85.4  9.12 162685124    30217  19.1         4.74   1501737720
 6 15.6   33.6 -86.8 10.2   26929603     9010  41.8        17.5    2878192209
 7 12.4   33.3 -87.0 10.2  166239542    16140   1.70        4.30   2878192209
 8 11.1   33.5 -87.3  8.16 385566685     3699   0           0.162  3423328940
 9 14.6   33.5 -86.9 10.2   10636977    11458  43.6        15.6    2878192209
10 12.0   33.7 -86.7  9.30 150661846    21725   1.48        4.25   2878192209
# … with 574 more rows, and 33 more variables: county_pop <dbl>,
#   log_dist_to_prisec <dbl>, log_pri_length_5000 <dbl>,
#   log_pri_length_25000 <dbl>, log_prisec_length_500 <dbl>,
#   log_prisec_length_1000 <dbl>, log_prisec_length_5000 <dbl>,
#   log_prisec_length_10000 <dbl>, log_prisec_length_25000 <dbl>,
#   log_nei_2008_pm10_sum_15000 <dbl>, log_nei_2008_pm10_sum_25000 <dbl>,
#   popdens_county <dbl>, popdens_zcta <dbl>, nohs <dbl>, somehs <dbl>,
#   hs <dbl>, somecollege <dbl>, associate <dbl>, bachelor <dbl>, grad <dbl>,
#   pov <dbl>, hs_orless <dbl>, urc2013 <dbl>, aod <dbl>,
#   state_Not.California <dbl>, city_Not.in.a.city <dbl>, .fitted <dbl>,
#   .se.fit <dbl>, .resid <dbl>, .hat <dbl>, .sigma <dbl>, .cooksd <dbl>,
#   .std.resid <dbl>

wflowoutput<-PM_wflow_fit %>% 
  pull_workflow_fit() %>% 
  broom::tidy() 


#The output is identical using workflows to fit the model or just parsnip
identical(tidy(PM_fit), wflowoutput)

[1] TRUE

OK, so we have fit our model on our training data, which means we have created a model to predict values of air pollution based on the predictors that we have included. Yay!

We can get a sense of the variable importance using the vip() function of the vip package.

Let’s take a look at the top 10 contributing variables:

PM_wflow_fit %>% 
  pull_workflow_fit() %>% 
  vip(num_features = 10)

Model Performance

Let’s take a look at how well our model fit our training data:

##using the parsnip version
#can grab the fitted values simply like this:
parsnip_fitted_values<-fitted(PM_fit[["fit"]])

parsnip_fitted_values

        1         2         3         4         5         6         7         8 
 9.461664 10.429189 11.795351 11.139746 10.863402 11.091857 10.073411  8.189021 
        9        10        11        12        13        14        15        16 
11.517888 10.618367  9.557499 10.218323  9.865608 11.011802 12.139040 10.202025 
       17        18        19        20        21        22        23        24 
10.517763  9.548571  6.837646  7.472453  7.018230 10.126795 10.133236  9.582305 
       25        26        27        28        29        30        31        32 
 9.702583  6.862089  7.005665  8.477788  8.892971  9.244236 11.417965 10.464960 
       33        34        35        36        37        38        39        40 
 9.652665  9.943594 11.435373 10.090206 10.304301 11.188076 11.193689 10.767763 
       41        42        43        44        45        46        47        48 
10.722817  9.949272 10.521367  9.976252 11.955651 11.674338 13.902418 13.582229 
       49        50        51        52        53        54        55        56 
14.713622 11.219505 11.320158 15.169192  8.682176 13.357477 11.166948 14.162085 
       57        58        59        60        61        62        63        64 
13.128209 10.867659 12.557864 11.404720 14.609938 10.003168 11.917995  9.886651 
       65        66        67        68        69        70        71        72 
11.161819 12.326977 12.722475 11.439546 11.380029 14.615337 10.945608 12.754387 
       73        74        75        76        77        78        79        80 
11.464256 13.050833 11.162822 14.092701 10.025276 12.181293 12.573405 12.850499 
       81        82        83        84        85        86        87        88 
12.523434 13.398889  8.867379 12.250365 15.385913 10.726820 12.320977 12.282267 
       89        90        91        92        93        94        95        96 
10.903510 13.504498 10.678776 11.367829 14.267153 11.959658 14.860756 12.245437 
       97        98        99       100       101       102       103       104 
11.009600 12.303933 10.376440  9.254959  8.458964 10.797185  9.615059  8.319766 
      105       106       107       108       109       110       111       112 
 6.584026  8.480982  7.889599  8.520672  9.507922  9.460726 12.291791 10.315367 
      113       114       115       116       117       118       119       120 
 9.899759 10.270368  9.531778  8.476645 10.506159 10.457589 11.548915 12.211470 
      121       122       123       124       125       126       127       128 
12.754005 11.760960 11.527420 13.358326 12.970527 11.450112  9.662891  8.860986 
      129       130       131       132       133       134       135       136 
 9.630399  9.107099  9.336214 10.377070 10.913731 10.012413 11.655663  9.197286 
      137       138       139       140       141       142       143       144 
 9.895577  8.992575 11.062182 10.928292 11.173204  9.290713  9.462250  9.815623 
      145       146       147       148       149       150       151       152 
10.567132 10.041853 12.438193 12.837043 11.155833 11.068409 11.568628 10.646408 
      153       154       155       156       157       158       159       160 
11.120305 10.674787 11.342914 11.600863  9.872202 12.933636 12.330424 10.674265 
      161       162       163       164       165       166       167       168 
10.692636  8.591527  8.483559  8.277465  9.219110  7.999004  7.325638 10.668789 
      169       170       171       172       173       174       175       176 
10.442171 13.278721 12.987068 12.609307 12.139812 13.311952  9.560588 11.231509 
      177       178       179       180       181       182       183       184 
 9.614572 11.388626 11.027188 10.701644 10.408321 11.281993 11.952935 12.314095 
      185       186       187       188       189       190       191       192 
11.343060 12.301626 11.670473 11.538465 11.450721 12.583606  9.908812 11.957886 
      193       194       195       196       197       198       199       200 
12.476480 11.356784 11.408551 11.386678 11.371265 11.869010 11.386166 10.716040 
      201       202       203       204       205       206       207       208 
13.547182 13.430563 14.488555 13.382721 11.721882 13.138956 11.617227 12.919960 
      209       210       211       212       213       214       215       216 
11.610622 11.234236 12.441517 11.565354 11.775386 11.870380 11.612343 10.513345 
      217       218       219       220       221       222       223       224 
11.332066 11.356102  8.622875  9.342198 10.532866 11.075768 10.814015  9.148910 
      225       226       227       228       229       230       231       232 
10.259892  8.851762 10.150785 11.702025 11.247041  9.448320  8.715848 10.465157 
      233       234       235       236       237       238       239       240 
 9.063738  8.512195  9.359375  8.694810 12.250888 10.950829 11.172992 11.113227 
      241       242       243       244       245       246       247       248 
12.539577 11.130235 11.520853 13.270299 12.445247 11.641586 11.084325 11.681937 
      249       250       251       252       253       254       255       256 
12.903454 11.675062 11.868743 11.887980 12.142425 11.125280 10.206552 13.083514 
      257       258       259       260       261       262       263       264 
12.553520 13.232212 11.590232 11.070216 11.072157 11.721994 10.595329 14.228343 
      265       266       267       268       269       270       271       272 
13.651194 10.111261 12.956021 13.017707 10.712920 13.277355 11.914351 11.578296 
      273       274       275       276       277       278       279       280 
12.601671 13.808564 13.137512 11.706520 11.860984  9.304048 10.849169 10.160520 
      281       282       283       284       285       286       287       288 
10.957044 10.853018 11.595010 11.208375 11.222850 10.646720 10.875621 10.875865 
      289       290       291       292       293       294       295       296 
12.139867  9.344793  7.829627  7.419011 10.276070  8.595271 11.984527 10.025066 
      297       298       299       300       301       302       303       304 
11.658014 10.844432 13.668283 13.877468 11.842378 12.609642 12.315624  7.550419 
      305       306       307       308       309       310       311       312 
10.543817 10.140751 10.640678  9.986984 10.672587 12.808831  8.608027  8.762374 
      313       314       315       316       317       318       319       320 
 9.800301  8.990358 10.913175 10.986842  9.848777 10.184492 10.730669  8.925857 
      321       322       323       324       325       326       327       328 
 9.810568 11.051639 12.301691 13.244583 12.143078  7.420562  7.181852  8.115063 
      329       330       331       332       333       334       335       336 
 7.557550  7.598764  7.963658  6.033631  8.174960  7.009403  6.547324  6.721532 
      337       338       339       340       341       342       343       344 
10.495751 10.127753  8.447220  8.955561  9.111061  6.229223  7.721007  8.702168 
      345       346       347       348       349       350       351       352 
 9.796073 10.735891 10.674564  9.993831  9.387035 11.404517 12.718761 12.490109 
      353       354       355       356       357       358       359       360 
12.490096 12.868844 14.159822 11.862943 12.490662 12.814318 13.292919 11.342429 
      361       362       363       364       365       366       367       368 
11.790078 13.543314 12.475061 11.439849  9.501710  9.938767  8.662747  9.361561 
      369       370       371       372       373       374       375       376 
 7.333849  7.619211  8.424590  8.182639  8.244907 11.421981 12.598988 12.766756 
      377       378       379       380       381       382       383       384 
10.004449 11.846089 11.484665  8.502844  9.985085 14.093239 11.022401 10.343070 
      385       386       387       388       389       390       391       392 
10.870999 11.580658  8.164514  9.816118 13.810603 10.765304 10.679352  8.943409 
      393       394       395       396       397       398       399       400 
11.394455 10.988983 11.911609 11.169336 10.926602 11.106484 10.975384 10.570735 
      401       402       403       404       405       406       407       408 
10.800890 11.216585 11.003865 12.058352  8.979768 10.417642 10.840787  8.564483 
      409       410       411       412       413       414       415       416 
10.468748  7.722598  9.024726  8.627277  8.712337 12.306801 11.259834 12.135047 
      417       418       419       420       421       422       423       424 
12.909887 13.900961 14.613681 14.659878 14.516719 11.834130 12.835139 11.524950 
      425       426       427       428       429       430       431       432 
11.205715 12.690535 12.366608 12.257143 11.554346 11.571765 10.268979 13.831592 
      433       434       435       436       437       438       439       440 
12.488541 12.274712 11.792411 12.159000 11.286455 12.398744 12.720250 11.839790 
      441       442       443       444       445       446       447       448 
12.042530 10.616712 11.592961  8.680015 10.002905 10.098331  9.303447  8.025186 
      449       450       451       452       453       454       455       456 
 8.160942  7.701500  7.573096  6.440509  8.413890  7.573817  8.654901  6.527707 
      457       458       459       460       461       462       463       464 
 9.316892  9.239795  7.577079  7.234717  7.744910 10.396131 11.349235 10.531142 
      465       466       467       468       469       470       471       472 
11.286622 11.324207 11.904075 11.984744 11.254344 10.087881 12.491983 12.727448 
      473       474       475       476       477       478       479       480 
11.023725  9.843285 12.532214 10.909423 11.556938 11.246725 13.678810 10.995575 
      481       482       483       484       485       486       487       488 
13.406587 12.947436 11.536939 10.519024 10.614975 10.691272 12.606993 11.966835 
      489       490       491       492       493       494       495       496 
11.220192  9.923407 11.012463  9.424555 10.077508 10.808538 10.687251 11.438914 
      497       498       499       500       501       502       503       504 
12.408781 10.434905  8.167398  7.968639  7.492743  7.971856  7.794390  8.631370 
      505       506       507       508       509       510       511       512 
 7.279254  7.833077 10.330354 10.925965 10.909906  8.503572 12.117403 10.667929 
      513       514       515       516       517       518       519       520 
10.129061 10.550327  9.918316 11.522068 11.296737 11.553654 10.476981 11.941762 
      521       522       523       524       525       526       527       528 
10.578826 11.080807  9.211676 10.323909  9.621231  8.335038  8.188955  9.082492 
      529       530       531       532       533       534       535       536 
 8.075345  9.420944  8.314572  9.247623  8.538591  9.859266  7.647124  9.902155 
      537       538       539       540       541       542       543       544 
 8.546524  9.920558 11.604894 11.918312 12.416602 10.819737 11.466657 11.179962 
      545       546       547       548       549       550       551       552 
11.665760  9.817793 14.172912 10.192052 11.355575 13.772991  8.847900  7.970775 
      553       554       555       556       557       558       559       560 
 7.460673  8.349559  8.135941 11.447727 11.618942 11.936942 10.918982 11.041702 
      561       562       563       564       565       566       567       568 
 9.949377 11.499303 10.904615 11.689807 11.747444 11.380585 11.009586 10.850464 
      569       570       571       572       573       574       575       576 
11.602371 11.551376 12.485626 11.270907 10.170153 10.619552  9.784335  8.033726 
      577       578       579       580       581       582       583       584 
10.810373  8.993775  5.743434  8.959191  8.089329  6.061699  8.339742  6.931294 

#can also get the values using the augment function
parsnip_fitted_values <- augment(PM_fit[["fit"]], data = juiced_train) %>% 
select(value, .fitted:.std.resid)

parsnip_fitted_values

# A tibble: 584 x 8
   value .fitted .se.fit .resid   .hat .sigma    .cooksd .std.resid
   <dbl>   <dbl>   <dbl>  <dbl>  <dbl>  <dbl>      <dbl>      <dbl>
 1  9.60    9.46   0.375  0.136 0.0324   2.09 0.00000433     0.0663
 2 10.8    10.4    0.383  0.371 0.0338   2.09 0.0000338      0.181 
 3 11.2    11.8    0.404 -0.583 0.0376   2.09 0.0000936     -0.285 
 4 12.4    11.1    0.388  1.24  0.0346   2.08 0.000384       0.604 
 5 10.5    10.9    0.426 -0.355 0.0418   2.09 0.0000388     -0.174 
 6 15.6    11.1    0.379  4.50  0.0332   2.08 0.00486        2.20  
 7 12.4    10.1    0.459  2.33  0.0485   2.08 0.00196        1.14  
 8 11.1     8.19   0.777  2.91  0.139    2.08 0.0108         1.51  
 9 14.6    11.5    0.378  3.03  0.0329   2.08 0.00219        1.48  
10 12.0    10.6    0.452  1.33  0.0470   2.08 0.000623       0.656 
# … with 574 more rows

##using the workflows version
wf_fit <-PM_wflow_fit %>% 
  pull_workflow_fit()

#can grab the fitted values simply like this:
wf_fitted_values<-fitted(wf_fit[["fit"]])

wf_fitted_values

        1         2         3         4         5         6         7         8 
 9.461664 10.429189 11.795351 11.139746 10.863402 11.091857 10.073411  8.189021 
        9        10        11        12        13        14        15        16 
11.517888 10.618367  9.557499 10.218323  9.865608 11.011802 12.139040 10.202025 
       17        18        19        20        21        22        23        24 
10.517763  9.548571  6.837646  7.472453  7.018230 10.126795 10.133236  9.582305 
       25        26        27        28        29        30        31        32 
 9.702583  6.862089  7.005665  8.477788  8.892971  9.244236 11.417965 10.464960 
       33        34        35        36        37        38        39        40 
 9.652665  9.943594 11.435373 10.090206 10.304301 11.188076 11.193689 10.767763 
       41        42        43        44        45        46        47        48 
10.722817  9.949272 10.521367  9.976252 11.955651 11.674338 13.902418 13.582229 
       49        50        51        52        53        54        55        56 
14.713622 11.219505 11.320158 15.169192  8.682176 13.357477 11.166948 14.162085 
       57        58        59        60        61        62        63        64 
13.128209 10.867659 12.557864 11.404720 14.609938 10.003168 11.917995  9.886651 
       65        66        67        68        69        70        71        72 
11.161819 12.326977 12.722475 11.439546 11.380029 14.615337 10.945608 12.754387 
       73        74        75        76        77        78        79        80 
11.464256 13.050833 11.162822 14.092701 10.025276 12.181293 12.573405 12.850499 
       81        82        83        84        85        86        87        88 
12.523434 13.398889  8.867379 12.250365 15.385913 10.726820 12.320977 12.282267 
       89        90        91        92        93        94        95        96 
10.903510 13.504498 10.678776 11.367829 14.267153 11.959658 14.860756 12.245437 
       97        98        99       100       101       102       103       104 
11.009600 12.303933 10.376440  9.254959  8.458964 10.797185  9.615059  8.319766 
      105       106       107       108       109       110       111       112 
 6.584026  8.480982  7.889599  8.520672  9.507922  9.460726 12.291791 10.315367 
      113       114       115       116       117       118       119       120 
 9.899759 10.270368  9.531778  8.476645 10.506159 10.457589 11.548915 12.211470 
      121       122       123       124       125       126       127       128 
12.754005 11.760960 11.527420 13.358326 12.970527 11.450112  9.662891  8.860986 
      129       130       131       132       133       134       135       136 
 9.630399  9.107099  9.336214 10.377070 10.913731 10.012413 11.655663  9.197286 
      137       138       139       140       141       142       143       144 
 9.895577  8.992575 11.062182 10.928292 11.173204  9.290713  9.462250  9.815623 
      145       146       147       148       149       150       151       152 
10.567132 10.041853 12.438193 12.837043 11.155833 11.068409 11.568628 10.646408 
      153       154       155       156       157       158       159       160 
11.120305 10.674787 11.342914 11.600863  9.872202 12.933636 12.330424 10.674265 
      161       162       163       164       165       166       167       168 
10.692636  8.591527  8.483559  8.277465  9.219110  7.999004  7.325638 10.668789 
      169       170       171       172       173       174       175       176 
10.442171 13.278721 12.987068 12.609307 12.139812 13.311952  9.560588 11.231509 
      177       178       179       180       181       182       183       184 
 9.614572 11.388626 11.027188 10.701644 10.408321 11.281993 11.952935 12.314095 
      185       186       187       188       189       190       191       192 
11.343060 12.301626 11.670473 11.538465 11.450721 12.583606  9.908812 11.957886 
      193       194       195       196       197       198       199       200 
12.476480 11.356784 11.408551 11.386678 11.371265 11.869010 11.386166 10.716040 
      201       202       203       204       205       206       207       208 
13.547182 13.430563 14.488555 13.382721 11.721882 13.138956 11.617227 12.919960 
      209       210       211       212       213       214       215       216 
11.610622 11.234236 12.441517 11.565354 11.775386 11.870380 11.612343 10.513345 
      217       218       219       220       221       222       223       224 
11.332066 11.356102  8.622875  9.342198 10.532866 11.075768 10.814015  9.148910 
      225       226       227       228       229       230       231       232 
10.259892  8.851762 10.150785 11.702025 11.247041  9.448320  8.715848 10.465157 
      233       234       235       236       237       238       239       240 
 9.063738  8.512195  9.359375  8.694810 12.250888 10.950829 11.172992 11.113227 
      241       242       243       244       245       246       247       248 
12.539577 11.130235 11.520853 13.270299 12.445247 11.641586 11.084325 11.681937 
      249       250       251       252       253       254       255       256 
12.903454 11.675062 11.868743 11.887980 12.142425 11.125280 10.206552 13.083514 
      257       258       259       260       261       262       263       264 
12.553520 13.232212 11.590232 11.070216 11.072157 11.721994 10.595329 14.228343 
      265       266       267       268       269       270       271       272 
13.651194 10.111261 12.956021 13.017707 10.712920 13.277355 11.914351 11.578296 
      273       274       275       276       277       278       279       280 
12.601671 13.808564 13.137512 11.706520 11.860984  9.304048 10.849169 10.160520 
      281       282       283       284       285       286       287       288 
10.957044 10.853018 11.595010 11.208375 11.222850 10.646720 10.875621 10.875865 
      289       290       291       292       293       294       295       296 
12.139867  9.344793  7.829627  7.419011 10.276070  8.595271 11.984527 10.025066 
      297       298       299       300       301       302       303       304 
11.658014 10.844432 13.668283 13.877468 11.842378 12.609642 12.315624  7.550419 
      305       306       307       308       309       310       311       312 
10.543817 10.140751 10.640678  9.986984 10.672587 12.808831  8.608027  8.762374 
      313       314       315       316       317       318       319       320 
 9.800301  8.990358 10.913175 10.986842  9.848777 10.184492 10.730669  8.925857 
      321       322       323       324       325       326       327       328 
 9.810568 11.051639 12.301691 13.244583 12.143078  7.420562  7.181852  8.115063 
      329       330       331       332       333       334       335       336 
 7.557550  7.598764  7.963658  6.033631  8.174960  7.009403  6.547324  6.721532 
      337       338       339       340       341       342       343       344 
10.495751 10.127753  8.447220  8.955561  9.111061  6.229223  7.721007  8.702168 
      345       346       347       348       349       350       351       352 
 9.796073 10.735891 10.674564  9.993831  9.387035 11.404517 12.718761 12.490109 
      353       354       355       356       357       358       359       360 
12.490096 12.868844 14.159822 11.862943 12.490662 12.814318 13.292919 11.342429 
      361       362       363       364       365       366       367       368 
11.790078 13.543314 12.475061 11.439849  9.501710  9.938767  8.662747  9.361561 
      369       370       371       372       373       374       375       376 
 7.333849  7.619211  8.424590  8.182639  8.244907 11.421981 12.598988 12.766756 
      377       378       379       380       381       382       383       384 
10.004449 11.846089 11.484665  8.502844  9.985085 14.093239 11.022401 10.343070 
      385       386       387       388       389       390       391       392 
10.870999 11.580658  8.164514  9.816118 13.810603 10.765304 10.679352  8.943409 
      393       394       395       396       397       398       399       400 
11.394455 10.988983 11.911609 11.169336 10.926602 11.106484 10.975384 10.570735 
      401       402       403       404       405       406       407       408 
10.800890 11.216585 11.003865 12.058352  8.979768 10.417642 10.840787  8.564483 
      409       410       411       412       413       414       415       416 
10.468748  7.722598  9.024726  8.627277  8.712337 12.306801 11.259834 12.135047 
      417       418       419       420       421       422       423       424 
12.909887 13.900961 14.613681 14.659878 14.516719 11.834130 12.835139 11.524950 
      425       426       427       428       429       430       431       432 
11.205715 12.690535 12.366608 12.257143 11.554346 11.571765 10.268979 13.831592 
      433       434       435       436       437       438       439       440 
12.488541 12.274712 11.792411 12.159000 11.286455 12.398744 12.720250 11.839790 
      441       442       443       444       445       446       447       448 
12.042530 10.616712 11.592961  8.680015 10.002905 10.098331  9.303447  8.025186 
      449       450       451       452       453       454       455       456 
 8.160942  7.701500  7.573096  6.440509  8.413890  7.573817  8.654901  6.527707 
      457       458       459       460       461       462       463       464 
 9.316892  9.239795  7.577079  7.234717  7.744910 10.396131 11.349235 10.531142 
      465       466       467       468       469       470       471       472 
11.286622 11.324207 11.904075 11.984744 11.254344 10.087881 12.491983 12.727448 
      473       474       475       476       477       478       479       480 
11.023725  9.843285 12.532214 10.909423 11.556938 11.246725 13.678810 10.995575 
      481       482       483       484       485       486       487       488 
13.406587 12.947436 11.536939 10.519024 10.614975 10.691272 12.606993 11.966835 
      489       490       491       492       493       494       495       496 
11.220192  9.923407 11.012463  9.424555 10.077508 10.808538 10.687251 11.438914 
      497       498       499       500       501       502       503       504 
12.408781 10.434905  8.167398  7.968639  7.492743  7.971856  7.794390  8.631370 
      505       506       507       508       509       510       511       512 
 7.279254  7.833077 10.330354 10.925965 10.909906  8.503572 12.117403 10.667929 
      513       514       515       516       517       518       519       520 
10.129061 10.550327  9.918316 11.522068 11.296737 11.553654 10.476981 11.941762 
      521       522       523       524       525       526       527       528 
10.578826 11.080807  9.211676 10.323909  9.621231  8.335038  8.188955  9.082492 
      529       530       531       532       533       534       535       536 
 8.075345  9.420944  8.314572  9.247623  8.538591  9.859266  7.647124  9.902155 
      537       538       539       540       541       542       543       544 
 8.546524  9.920558 11.604894 11.918312 12.416602 10.819737 11.466657 11.179962 
      545       546       547       548       549       550       551       552 
11.665760  9.817793 14.172912 10.192052 11.355575 13.772991  8.847900  7.970775 
      553       554       555       556       557       558       559       560 
 7.460673  8.349559  8.135941 11.447727 11.618942 11.936942 10.918982 11.041702 
      561       562       563       564       565       566       567       568 
 9.949377 11.499303 10.904615 11.689807 11.747444 11.380585 11.009586 10.850464 
      569       570       571       572       573       574       575       576 
11.602371 11.551376 12.485626 11.270907 10.170153 10.619552  9.784335  8.033726 
      577       578       579       580       581       582       583       584 
10.810373  8.993775  5.743434  8.959191  8.089329  6.061699  8.339742  6.931294 

#can also get the values using the augment function
wf_fitted_values <-augment(wf_fit[["fit"]], data = juiced_train) %>% 
select(value, .fitted:.std.resid)

wf_fitted_values

# A tibble: 584 x 8
   value .fitted .se.fit .resid   .hat .sigma    .cooksd .std.resid
   <dbl>   <dbl>   <dbl>  <dbl>  <dbl>  <dbl>      <dbl>      <dbl>
 1  9.60    9.46   0.375  0.136 0.0324   2.09 0.00000433     0.0663
 2 10.8    10.4    0.383  0.371 0.0338   2.09 0.0000338      0.181 
 3 11.2    11.8    0.404 -0.583 0.0376   2.09 0.0000936     -0.285 
 4 12.4    11.1    0.388  1.24  0.0346   2.08 0.000384       0.604 
 5 10.5    10.9    0.426 -0.355 0.0418   2.09 0.0000388     -0.174 
 6 15.6    11.1    0.379  4.50  0.0332   2.08 0.00486        2.20  
 7 12.4    10.1    0.459  2.33  0.0485   2.08 0.00196        1.14  
 8 11.1     8.19   0.777  2.91  0.139    2.08 0.0108         1.51  
 9 14.6    11.5    0.378  3.03  0.0329   2.08 0.00219        1.48  
10 12.0    10.6    0.452  1.33  0.0470   2.08 0.000623       0.656 
# … with 574 more rows

#gives us the same fitted values
identical(parsnip_fitted_values, wf_fitted_values)

[1] TRUE

## Let's make a plot of fitted and real values

ggplot(wf_fitted_values, aes(x = .fitted, y = value)) + geom_point()

Ok, so our fitted range appears to be smaller than the real values. We could probably do a bit better.

Let’s take a look at how well our model seems to be preforming more formally:

When assessing the performance of a model, the metrics we use depend on if we are preforming a classifcation or prediction also known as regression analysis. In our case we are performing a regression or prediction analysis and the metrics often used are: 1) mean absolute error (mae)
2) R squared error (rsq) This is also known as the coefficient of determiniation which is the squared correlation between truth and estimate
3) root mean squared error (rmse)

We can use the yardstick package to quickly calculate estimates for all of these values using the metrics() function. Alternatively if you only wanted one metric you could use the mae(), rsq(), or rmse() functions respectively. This is helpful to examine with our fitted training set values to see how well our model is performing and if we need to make adjustments.

yardstick::metrics(wf_fitted_values, truth = value, estimate = .fitted)

# A tibble: 3 x 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard       2.02 
2 rsq     standard       0.433
3 mae     standard       1.49 

yardstick::mae(wf_fitted_values, truth = value, estimate = .fitted)

# A tibble: 1 x 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 mae     standard        1.49

Cross validation sample splitting

We will use the rsample package again in order to further implement what are called cross validation techniques. This is also called resampling or repartioning. Note: we are not actually getting new samples from the underlying distribution so the term resampling is a bit of a misnomer.

Cross validation splits our training data into multiple training data sets to allow for a deeper assessment of the accuracy of the model.

Here is a visualization of the concept for cross validation/resampling/repartitioning from Max Kuhn:

Technically creating our testing and training set out of our original training data is sometimes considered a form of cross validation called the holdout method. As we just learned this can give us a better sense of the accuracy of our data in a more generalizable way.

However, we can do a better job of optimizing our model for accuracy if we also perform another type of cross validation on the newly defined training set that we just created. There are many cross validation methods and most can be easily implemented using rsamples package. We will use a very popular method called either k-fold or v-fold cross validation.

This method involves essentially preforming the hold out method iteratively with the training data.

First the training set is divided into k or v equally sized smaller pieces.

Then the model is trained on the model on k-1 or v-1 subsets of the data iteratively (removing a different v or k untill all possible k-1 or v-1 sets have been evaluated) to get a sense of the performance of the model. This is really useful for fine tuning specific aspects of the model in a process called model tuning.

Here is a visualization of how the folds are created:

Note: People typically ignore spatial dependence with cross validation of air pollution monitoring data in the air pollution field, so we will do the same. However, it might make sense to leave out blocks of monitors rather than random individual monitors to help account for some spatial dependence.

The vfold_cv() function of the rsample package can be used to parse the training data into folds for k-fold/v-fold cross validation.

The v argument specifies the number of folds to create. The repeats argument specifies if any samples should be repeated across folds - defualt is FALSE The strata argument specifies a variable to stratify samples across folds (just like in initial_split()).

Again because these are created at random, we need to use the base set.seed() function in order to obtain the same results each time we knit this document. Generally speaking using 10 folds is good practice, but this depends on the variablity within your data. We are going to use 4 for the sake of expediency.

set.seed(1234)

vfold_pm <-rsample::vfold_cv(data = train_pm, v = 4)

vfold_pm

#  4-fold cross-validation 
# A tibble: 4 x 2
  splits            id   
  <list>            <chr>
1 <split [438/146]> Fold1
2 <split [438/146]> Fold2
3 <split [438/146]> Fold3
4 <split [438/146]> Fold4

vfold_pm$splits$`1`

NULL

vfold_pm$splits$`2`

NULL

Once the folds are created they can be used to evaluate peformance by fitting the model to each of the resamples that we created:

We can fit the model to our cross validation folds using the fit_resamples() function of the tune package, by specifying our workflow object and the cross validation fold object we just created. See here for more information.

control <- control_resamples(save_pred = TRUE)
resample_fit <-tune::fit_resamples(PM_wflow, vfold_pm)

We can now take a look at various metrics of performance based on the fit of our cross validation “resamples”. To do this we will use the show_best() function of the tune package.

tune::show_best(resample_fit, metric = "rmse", "rsq")

# A tibble: 1 x 5
  .metric .estimator  mean     n std_err
  <chr>   <chr>      <dbl> <int>   <dbl>
1 rmse    standard    2.18     4  0.0763

Tuning

Now let’s try some tuning.

Let’s take a closer look at how the air pollution monitor values vary with the location latitude and longitude.

train_pm %>% 
  dplyr::select(value, lon, lat) %>% 
  tidyr::pivot_longer(cols = c(lon, lat), 
                      names_to = "predictor", values_to = "loc_value") %>% 
  ggplot(aes(x = loc_value, value)) + 
  geom_point(alpha = .2) + 
  geom_smooth(se = FALSE) + 
  #scale_y_log10() +
  facet_wrap(~ predictor, scales = "free_x")

We can see that there does not appear to be a single linear relationship for either of these predictors. Thus we might want to think about using splines or this(https://towardsdatascience.com/numerical-interpolation-natural-cubic-spline-52c1157b98ac) or just this(https://tidymodels.github.io/tune/articles/getting_started.html) or this(https://www.psych.mcgill.ca/misc/fda/ex-basis-b1.html) to model the relationship in our training data more closely. For example for the latitude plot (left) if we had 2 lines and one breakpoint called a knot around 40, with the first line having a positive slope and the second with a negative slope this would fit the data more similarly to the blue line shown in the figure.

We can tune for the number of knots by using a step function in the recipes package called step_ns() where ns stands for natural splines. In order to tune for the number of knots or degrees of freedom, we can set the deg_free argument to tune(). This is helpful, becuase we aren’t exactly sure how closely we should be following the relationship with the value and our longitude and latitude data in our training data to achieve good accuracy yet keep our model generalizable for other data.

This is when our cross validation methods become really handy. We can test out different values for the deg_free argument and see how our model performance varies across our training folds to try to find the optimal value.

We will update our recipe to add these steps. It is a good idea to do this for individual predictors because you can name each with the tune argument so that you can keep track of it later. We can see what we intend to tune with the parameters() function of the dials package.

See here for more information about implementing this in tidymodels.

novel_rec %<>%
  step_ns(lon, deg_free = tune("lon df")) %>%
  step_ns(lat,  deg_free = tune("lat df"))
# novel_rec %<>%
#   step_ns(lat,  deg_free = tune())
# novel_rec

pm_param <-dials::parameters(novel_rec)
pm_param

Collection of 2 parameters for tuning

     id parameter type object class
 lon df       deg_free    nparam[+]
 lat df       deg_free    nparam[+]

Genearlly you could use the grid_*() functions of the dials package to create the different combinations of degrees of freedom to test for both variables to optimize the model. In our case we can visably see that if we add more than say 4 or 5 degrees of freedom we will likely overfit the data. So instead of using these functions we will create our own grid using the base seq() and expand.grid() functions.

#an example of what you could do:
#spline_grid <-dials::grid_regular(pm_param, levels = 3)
df_vals <- seq(1, 5, by = 2)
spline_grid <- expand.grid(`lon df` = df_vals, `lat df` = df_vals)
spline_grid

  lon df lat df
1      1      1
2      3      1
3      5      1
4      1      3
5      3      3
6      5      3
7      1      5
8      3      5
9      5      5

Now we will tune this hyper-parameter (degrees of freedom) for both the lat and lon variables using our cross validation folds. To do this we will use the tune_grid() function of the tune package.

df_tuning <-lm_PM_model %>% 
  tune::tune_grid(novel_rec, resamples =vfold_pm, 
                  grid = spline_grid)

#df_tuning <-PM_wflow %>% tune::tune_grid(resamples =vfold_pm, 
 #                                        grid = spline_grid, 
  
#                                       param_set =pm_param)

df_tuning

#  4-fold cross-validation 
# A tibble: 4 x 4
  splits            id    .metrics          .notes           
  <list>            <chr> <list>            <list>           
1 <split [438/146]> Fold1 <tibble [18 × 5]> <tibble [18 × 1]>
2 <split [438/146]> Fold2 <tibble [18 × 5]> <tibble [18 × 1]>
3 <split [438/146]> Fold3 <tibble [18 × 5]> <tibble [18 × 1]>
4 <split [438/146]> Fold4 <tibble [18 × 5]> <tibble [18 × 1]>

df_tuning %>%
  collect_metrics()

# A tibble: 18 x 7
   `lon df` `lat df` .metric .estimator  mean     n std_err
      <dbl>    <dbl> <chr>   <chr>      <dbl> <int>   <dbl>
 1        1        1 rmse    standard   2.18      4  0.0763
 2        1        1 rsq     standard   0.361     4  0.0486
 3        1        3 rmse    standard   2.13      4  0.0807
 4        1        3 rsq     standard   0.383     4  0.0388
 5        1        5 rmse    standard   2.13      4  0.0812
 6        1        5 rsq     standard   0.386     4  0.0403
 7        3        1 rmse    standard   2.10      4  0.0613
 8        3        1 rsq     standard   0.402     4  0.0296
 9        3        3 rmse    standard   2.05      4  0.0709
10        3        3 rsq     standard   0.428     4  0.0191
11        3        5 rmse    standard   2.03      4  0.0656
12        3        5 rsq     standard   0.439     4  0.0197
13        5        1 rmse    standard   2.11      4  0.0616
14        5        1 rsq     standard   0.397     4  0.0258
15        5        3 rmse    standard   2.05      4  0.0739
16        5        3 rsq     standard   0.427     4  0.0164
17        5        5 rmse    standard   2.02      4  0.0717
18        5        5 rsq     standard   0.442     4  0.0153

show_best(df_tuning, metric = "rmse", n =1)

# A tibble: 1 x 7
  `lon df` `lat df` .metric .estimator  mean     n std_err
     <dbl>    <dbl> <chr>   <chr>      <dbl> <int>   <dbl>
1        5        5 rmse    standard    2.02     4  0.0717

Linear Regression Model with PCA

We can create another workflow to see how model pefromance compares using a different model. In this case we are going to perform something called Principal Component Analysis or PCA.

So what is PCA?

PCA is a widely used dimensionality reduction method (a form of unspervised machine learning). It creates new variables that capture the most variation within the data, yet reduce the data down to just a number of principal components. It does so by transforiming the data using orthogonal%20%3D%200.){target="_blank"} linear transformation. In otherwords, it creates new variables that are linear combinations of the variables within the data. Importantly these new variables are orthogonal, meaning that the new variables have zero covariance. In simpler terms, we are expessing unique types of variation within the data as new variables.

Check out this video for more information.

lm_PCA_rec <-recipe(train_pm) %>%
    update_role(everything(), new_role = "predictor")%>%
    update_role(value, new_role = "outcome")%>%
    update_role(id, new_role = "id variable") %>%
    update_role("fips", new_role = "county id") %>%
    step_dummy(state, county, city, zcta) %>%
    step_pca(all_predictors()) 

Let’s take a look to see what the step_pca function does to our predictors. To do so recall that we need to use the prep and juice functions of the recipes package on our recipe.

prepped_rec <- prep(lm_PCA_rec, verbose = TRUE, retain = TRUE )

oper 1 step dummy [training] 
oper 2 step pca [training] 
The retained training set is ~ 0.13 Mb  in memory.

juiced_train<- juice(prepped_rec)
glimpse(juiced_train)

Rows: 584
Columns: 8
$ id    <fct> 1003.001, 1027.0001, 1033.1002, 1055.001, 1069.0003, 1073.0023,…
$ value <dbl> 9.597647, 10.800000, 11.212174, 12.375394, 10.508850, 15.591017…
$ fips  <fct> 1003, 1027, 1033, 1055, 1069, 1073, 1073, 1073, 1073, 1073, 107…
$ PC1   <dbl> -4120232815, -1570558844, -1534934869, -1388102571, -1504354268…
$ PC2   <dbl> -118909004, -346706197, 10140574, -129802253, -136385028, 23433…
$ PC3   <dbl> -72827.994, 7132.812, -56957.316, 45519.861, 36487.258, 448189.…
$ PC4   <dbl> -20385.63273, 1615.38773, -7187.37037, -15876.58640, -25852.952…
$ PC5   <dbl> -846.027447, 90.287809, 171.426032, -576.131718, -926.823714, -…

We still want to use the lm package for our regression so we can use the same model object as before:

lm_PM_model

Linear Regression Model Specification (regression)

Computational engine: lm 

pca_wflow <-workflows::workflow() %>%
            workflows::add_recipe(lm_PCA_rec) %>%
            workflows::add_model(lm_PM_model)


pca_wflow

══ Workflow ═══════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: linear_reg()

── Preprocessor ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 Recipe Steps

● step_dummy()
● step_pca()

── Model ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Linear Regression Model Specification (regression)

Computational engine: lm 

Remember that using workflows we dont actually need to prep our recipe, we can just fit our model directly.

Fit the cross validation samples:

resample_pca_fit <-tune::fit_resamples(pca_wflow, vfold_pm)

Look at the performance:

collect_metrics(resample_pca_fit)

# A tibble: 2 x 5
  .metric .estimator   mean     n std_err
  <chr>   <chr>       <dbl> <int>   <dbl>
1 rmse    standard   2.76       4  0.271 
2 rsq     standard   0.0527     4  0.0184

And we can compare this with our previous performance:

collect_metrics(resample_fit)

# A tibble: 2 x 5
  .metric .estimator  mean     n std_err
  <chr>   <chr>      <dbl> <int>   <dbl>
1 rmse    standard   2.18      4  0.0763
2 rsq     standard   0.361     4  0.0486

So we can see that our performance isn’t quite as good - especially if we look at the rmse value.

Random Forest

Now for one last recipe, we are going to predict using a decision tree method called random forests.

A decision tree is a tool to partition data based on a series of sequential binary decisions.

knitr::include_graphics("https://miro.medium.com/max/1000/1*LMoJmXCsQlciGTEyoSN39g.jpeg")

[source]

knitr::include_graphics("https://miro.medium.com/max/1400/0*f_qQPFpdofWGLQqc.png")

In our case, the random forest algorithm that we are working with does not work well when there are categorical variabels with more than 53 levels, so we will need to remove the zcta variable.

RF_rec <-recipe(train_pm) %>%
    update_role(everything(), new_role = "predictor")%>%
    update_role(value, new_role = "outcome")%>%
    update_role(id, new_role = "id variable") %>%
    update_role("fips", new_role = "county id") %>%
    step_string2factor("state", "county", "city") %>%
    step_rm("county") %>%
    step_rm("id") %>%
    step_rm("fips")%>%
    step_rm("zcta") %>%
    step_corr(all_numeric())%>%
    step_nzv(all_numeric())

PMtree_model <- parsnip::rand_forest(mtry = 10, min_n = 3, mode = "regression")
PMtree_model

Random Forest Model Specification (regression)

Main Arguments:
  mtry = 10
  min_n = 3

RF_PM_model <- 
  PMtree_model  %>%
  parsnip::set_engine("randomForest") %>%
  set_mode("regression")

RF_PM_model

Random Forest Model Specification (regression)

Main Arguments:
  mtry = 10
  min_n = 3

Computational engine: randomForest 

RF_wflow <-workflows::workflow() %>%
            workflows::add_recipe(RF_rec) %>%
            workflows::add_model(RF_PM_model)
RF_wflow

══ Workflow ═══════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: rand_forest()

── Preprocessor ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
7 Recipe Steps

● step_string2factor()
● step_rm()
● step_rm()
● step_rm()
● step_rm()
● step_corr()
● step_nzv()

── Model ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Random Forest Model Specification (regression)

Main Arguments:
  mtry = 10
  min_n = 3

Computational engine: randomForest 

Fitting the data with just parsnip and with the workflow:

##just parsnip
RF_fit <- RF_PM_model %>% 
parsnip::fit(value ~., data =juiced_train_ready)

## with workflow
RF_wflow_fit <- parsnip::fit(RF_wflow, data = train_pm)
RF_wflow_fit

══ Workflow [trained] ═════════════════════════════════════════════════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: rand_forest()

── Preprocessor ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
7 Recipe Steps

● step_string2factor()
● step_rm()
● step_rm()
● step_rm()
● step_rm()
● step_corr()
● step_nzv()

── Model ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

Call:
 randomForest(x = as.data.frame(x), y = y, mtry = ~10, nodesize = ~3) 
               Type of random forest: regression
                     Number of trees: 500
No. of variables tried at each split: 10

          Mean of squared residuals: 3.016711
                    % Var explained: 58.14

Let’s take a look at the top 10 contributing variables:

RF_wflow_fit%>% 
  pull_workflow_fit() %>% 
  vip(num_features = 10)

Interesting, in the previous model the CMAQ values were also important, however the variable about if the monitor was located in California or not was also very predictive.

Now let’s take a look at model performance by fitting the data using cross validation:

resample_RF_fit <-tune::fit_resamples(RF_wflow, vfold_pm)
collect_metrics(resample_RF_fit)

# A tibble: 2 x 5
  .metric .estimator  mean     n std_err
  <chr>   <chr>      <dbl> <int>   <dbl>
1 rmse    standard   1.79      4  0.105 
2 rsq     standard   0.580     4  0.0377

Now let’s compare the performance of this model with the others:

# our initial linear regression model:
collect_metrics(resample_fit)

# A tibble: 2 x 5
  .metric .estimator  mean     n std_err
  <chr>   <chr>      <dbl> <int>   <dbl>
1 rmse    standard   2.18      4  0.0763
2 rsq     standard   0.361     4  0.0486

# our initial linear regression model with lat/lon degrees of freedom tuning:
show_best(df_tuning, metric = "rmse", n =1)

# A tibble: 1 x 7
  `lon df` `lat df` .metric .estimator  mean     n std_err
     <dbl>    <dbl> <chr>   <chr>      <dbl> <int>   <dbl>
1        5        5 rmse    standard    2.02     4  0.0717

# our PCA model:
collect_metrics(resample_pca_fit)

# A tibble: 2 x 5
  .metric .estimator   mean     n std_err
  <chr>   <chr>       <dbl> <int>   <dbl>
1 rmse    standard   2.76       4  0.271 
2 rsq     standard   0.0527     4  0.0184

Ok, so our first model had a mean rmse value of 2.18. The model with the lat/long degrees of freedom tuning had a mean rmse value of 2.02, thus showing some improvement. The PCA model had a mean rmse value of 2.76.

It looks like the random forest model had the lowest rmse value of 1.79.

If we tuned our random forest model based on the number of trees or the value for mtry (which is “The number of predictors that will be randomly sampled at each split when creating the tree models”), we might get a model with even better performance.

However, our cross valided mean rmse value of 1.79 is quite good because our range of true outcome values is much larger:3.4963303, 22.259123.

Final Model Performance Evaluation

Now that we have decided that we have resonable performance with the training data, we could stop here and use the yardstick package (and tune if using workflows to fit our model) to evaluate performance with our testing data.

So now we will use our random forests model to predict values for the monitors in the testing data.

Using parsnip we would need to use the baked data testing data. With the workflows package, we could use the raw testing data.

Importantly, ID variables are not delt with as nicely as with the workflows package so we would need to remove them. We did this above when created the processed training data for this model, the juiced_train data as well.

 baked_test_pm_ready <-baked_test_pm %>%select( -"id", -"fips")
 values_pred_parsnip <-predict(RF_fit, baked_test_pm_ready)
 values_pred_parsnip

# A tibble: 292 x 1
   .pred
   <dbl>
 1 11.1 
 2 11.8 
 3 12.0 
 4 11.3 
 5 11.8 
 6 11.9 
 7 10.6 
 8 10.6 
 9  8.65
10  7.60
# … with 282 more rows

# using the workflows version
values_pred_wfs <- 
  predict(RF_wflow_fit, test_pm) %>% 
  bind_cols(test_pm %>% select(value, fips, county, id)) 
values_pred_wfs

# A tibble: 292 x 5
   .pred value fips  county     id       
   <dbl> <dbl> <fct> <chr>      <fct>    
 1 11.1  11.7  1049  DeKalb     1049.1003
 2 12.0  13.1  1073  Jefferson  1073.101 
 3 11.6  12.2  1073  Jefferson  1073.2006
 4 11.3  12.2  1089  Madison    1089.0014
 5 11.5  11.4  1103  Morgan     1103.0011
 6 12.0  12.2  1121  Talladega  1121.0002
 7 10.8  10.9  4013  Maricopa   4013.4003
 8 10.7  10.6  4021  Pinal      4021.0001
 9  8.60 14.1  4023  Santa Cruz 4023.0004
10  7.59  5.83 4025  Yavapai    4025.2002
# … with 282 more rows

#model performance with test set
yardstick::metrics(values_pred_wfs, truth = value, estimate = .pred)

# A tibble: 3 x 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard       1.49 
2 rsq     standard       0.612
3 mae     standard       1.06 

Awesome! We can see that our rmse of 1.49 is quite similar with our testing data. We achieved quite good performance, which suggests that we would could predict other locations with more sparse monitoring based on our predictors with resonable accuracy.

We could also use the last_fit() function of the tune package to look at performance if we chose to create a workflow using the workflows package.

overallfit <-tune::last_fit(RF_wflow, pm_split)
 # or
overallfit <-RF_wflow %>%
  tune::last_fit(pm_split)

overallfit$.metrics[[1]]

# A tibble: 2 x 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard       1.49 
2 rsq     standard       0.607

We could check out test performance using the collect_metrics() function of the tune package.

test_performance <- overallfit  %>%tune::collect_metrics()
test_performance

# A tibble: 2 x 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard       1.49 
2 rsq     standard       0.607

Here you can see the predictions for the test set (the 292 rows with predicitions out of the 876 original monitor values) also using the tune package with the collect_predictions() function.

test_predictions <- overallfit  %>%tune::collect_predictions()
test_predictions

# A tibble: 292 x 4
   id               .pred  .row value
   <chr>            <dbl> <int> <dbl>
 1 train/test split 11.1      4 11.7 
 2 train/test split 11.9     10 13.1 
 3 train/test split 11.8     12 12.2 
 4 train/test split 11.4     15 12.2 
 5 train/test split 11.5     19 11.4 
 6 train/test split 12.1     22 12.2 
 7 train/test split 10.6     30 10.9 
 8 train/test split 10.7     36 10.6 
 9 train/test split  8.54    39 14.1 
10 train/test split  7.46    40  5.83
# … with 282 more rows

Data Visualization

Our main question for this case study was:

1. Can we predict annual average air pollution concentrations at the granularity of zip code regional levels using predictors such as data about population density, urbanization, road density, as well as, satellite pollution data and chemical modeling data?

We have indeed created a model that can predict fine particulate matter air pollution levels based on our predictor variables.

Now let’s make a plot of our predicted values and the true values.

First, let’s start by making a plot of our monitors:

We will use the following packages to create a map of the US: 1)sf 2)rnaturalearth 3)rnaturalearthdata 4)rgeos

According to this link on wikipedia, these are the latitude and longitude bounds of the contintental US.

top = 49.3457868 # north lat left = -124.7844079 # west long right = -66.9513812 # east long bottom = 24.7433195 # south lat

We will start with getting an outline of the US with the ne_countries() function of the rnaturalearth package.

library("sf")
library("rnaturalearth")
#library("rnaturalearthdata")
library("rgeos")

world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)

[1] "sf"         "data.frame"

ggplot(data = world) +
    geom_sf() +
    geom_point(data = pm, aes(x = lon, y = lat), size = 2, 
        shape = 23, fill = "darkred") +
    coord_sf(xlim = c(-125, -66), ylim = c(24.5, 50), expand = FALSE)

Now let’s add county lines.

County data is available from the maps package.

counties <- st_as_sf(maps::map("county", plot = FALSE, fill = TRUE))

monitors <-ggplot(data = world) +
   geom_sf(data = counties, fill = NA, color = gray(.5))+
    geom_point(data = pm, aes(x = lon, y = lat), size = 2, 
        shape = 23, fill = "darkred") +
    coord_sf(xlim = c(-125, -66), ylim = c(24.5, 50), expand = FALSE)+
  ggtitle("Monitor Locations")

monitors

Now let’s add a fill at the county level for the true monitor values of air pollution:

library(lwgeom)
library(stringr)

head(counties)

Simple feature collection with 6 features and 1 field
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -88.01778 ymin: 30.24071 xmax: -85.06131 ymax: 34.2686
CRS:            EPSG:4326
               ID                           geom
1 alabama,autauga MULTIPOLYGON (((-86.50517 3...
2 alabama,baldwin MULTIPOLYGON (((-87.93757 3...
3 alabama,barbour MULTIPOLYGON (((-85.42801 3...
4    alabama,bibb MULTIPOLYGON (((-87.02083 3...
5  alabama,blount MULTIPOLYGON (((-86.9578 33...
6 alabama,bullock MULTIPOLYGON (((-85.66866 3...

pm <-readr::read_csv(here("docs", "pm25_data.csv"))
counties %<>% tidyr::separate(ID, into = c("state", "county"), sep = ",")

counties[["county"]] <-str_to_title(counties[["county"]])

map_data <-inner_join(counties, pm, by = "county")

truth<-ggplot(data = world) +
    geom_sf() +
    geom_sf(data = map_data, aes(fill = value)) +
    #scale_fill_viridis_c(trans = "sqrt", alpha = .4) +
      #geom_point(data = pm, aes(x = lon, y = lat), size = 1.5, alpha = 2,
       # shape = 23, fill = "darkred") +
        coord_sf(xlim = c(-125, -66), ylim = c(24.5, 50), expand = FALSE)+
      ggtitle("Monitor PM2.5 levels")
truth

Now let’s do the same with our predicted values.

Let’s grab both the testing and training fitted values so that we have as much data as possible. In this case, the output structure for the training data fit is slightly different using randomForest. The fitted values are called predicted and the broom functions like tidy() and augment() will not work. So we will manually grab the fitted training data values.

#test data
values_pred_wfs

# A tibble: 292 x 5
   .pred value fips  county     id       
   <dbl> <dbl> <fct> <chr>      <fct>    
 1 11.1  11.7  1049  DeKalb     1049.1003
 2 12.0  13.1  1073  Jefferson  1073.101 
 3 11.6  12.2  1073  Jefferson  1073.2006
 4 11.3  12.2  1089  Madison    1089.0014
 5 11.5  11.4  1103  Morgan     1103.0011
 6 12.0  12.2  1121  Talladega  1121.0002
 7 10.8  10.9  4013  Maricopa   4013.4003
 8 10.7  10.6  4021  Pinal      4021.0001
 9  8.60 14.1  4023  Santa Cruz 4023.0004
10  7.59  5.83 4025  Yavapai    4025.2002
# … with 282 more rows

#training data
training_RF_estimates <-RF_fit[["fit"]][["predicted"]]
training_RF_estimates <-tibble(".pred" =training_RF_estimates,
value= train_pm$value, fips = train_pm$fips,  county = train_pm$county, id =train_pm$id)

all_pred <-bind_rows(values_pred_wfs, training_RF_estimates )

map_data <-inner_join(counties, all_pred, by = "county")

pred <-ggplot(data = world) +
    geom_sf() +
    geom_sf(data = map_data, aes(fill = .pred)) +
    #scale_fill_viridis_c(trans = "sqrt", alpha = .4) +
      #geom_point(data = pm, aes(x = lon, y = lat), size = 1.5, alpha =2,
      #  shape = 23, fill = "darkred") +
        coord_sf(xlim = c(-125, -66), ylim = c(24.5, 50), expand = FALSE)+
  ggtitle("Predicted PM2.5 levels")

pred

library(cowplot)

cowplot::plot_grid(monitors, truth, pred)

Summary

Let’s review everything:

We have explored gravametric fine particulate matter air pollution monitoring data. We have utilized the tidymodels ecosystem of packages to predict monitor values using a variety of predictors, also known as explanatory varaibles, including sattelite data, road density data, and population density, among others. Our model could be extended to be used to predict pollution levels in areas without monitoring, to help identify regions where populations maybe especially at risk for the health effects of air pollution.

We learned that there are two major types of what is called supervised machine learning: prediction and classification. We learned that prediction is used when the outcome variable is numeric and classifaction is performed when the outcome variable is categorical.

We performed the major steps of machine learning that we introduced in the beginning of the data analysis:

1. Data exploration

We used a packages like skimr, summarytools, corrplot, ggcorrplot, and GGally to better understand our data. These packages gave can tell us how many missing values each variable has (if any), the class of each variable, the distribution of values for each variable, the sparsity of each variable, and the level of correlation between variables.

2. Data splitting

We used the rsample package to first perform an initial split of our data into two peices: a training set and a testing set. The training set was used to optimize the model, while the testing set was used only to evaluate the performance of our final model. We also used the rsample package to create cross validation subsets of our training data. This allowed us to better assess the performance of our tested models using our training data.

3. Variable assignment and pre-processing

We used the recipes package to assign variable roles (such as outcome, predictor, and id variable). We also used this package to create a recipe for pre-processing our training and testing data. This involved steps such as: step_dummy to create dummy numeric encodings of our categorical variables, step_corr to remove highly correlated variables, step_nzv to remove near zero variance variabels that would contibute little to our model and potentially add noise. We learned that once our recipe was created and prepped using prep()we could extract the pre-processed training data using juice() or our pre-processed testing data using bake(). We also learned that if we used the newer workflows package that we did not need to the prep(), juice(), or bake() functions, but that it is still useful to know how to do so if we want to look at our data and how the recipe is influencing it more deeply.

4. Model specification, fitting, tuning and performance evaluation using the training data

We learned that the model needs to first be fit to the training data. We learned that in both classification and prediction, the model is fit to the traing data and the explanatory variables are used to estimate numeric values (in the case of prediction) or categorical values (in the case of classification) of the outcome variable of interest. We learned that we specify the model and its specifications using the parnsip package and that we also use this packae to fit the model using the fit() function. We learned that we if just use parsnip to fit the model, then we need to use the pre-processed training data (output from juice()). We learned that we can use the raw training data if we use the workflows package to create a workflow that pre-processes our data for us.

We learned that if the model fits well than the estimated values will be very similar to the true outcome variable values in our training data. We learned that we can assess model peformance using the yardstick package and the metrics() function. We also learned that we can use subsets of our training data (which we created with the rsample package) to perform cross validation to get a better estimate about the performance of our model using our training data, as we want our results to be generalizable and to perform well with other data, not just our training data. We used the fit_resamples() function of the tune package to fit our model on our different training data subsets and the collect_metrics() function (also of the tune package) to evaluate model performance using these subsets. We also learned that we can potentially improve model performance by tuning aspects about the model called hyperparameters to determine the best option for model performance. We learned that we can do this using the tune and dials packages and evaluating the performance of our model with the different hyperparameter options and our training data subsets that we used for cross validation. After we tested several different methods to model our data, we compared them to choose the best performing model as our final model.

5. Overall model performance evaluation

Once we chose our final model, we evaluated the final model performance using the testing data. This gives us a better estimate about how well the model will predict or classify the outcome variable of interest with new independent data. Ideally one would also perform an evaluation with independent data to provide a sense of how generalizable the model is to other data sources.

We first fit our model to our testing data using either just parsnip and the pre-processed testing data (using the bake() recipes function), or our raw testing data if we used a workflow. We used the same performance evaluation functions (yardstick::metrics() and tune::collect_metrics()(when using cross validion)). We also learned how we can use the last_fit() function of the tune package if we used a workflow to get the test data performance using the initial data and the testing/training split information.

Analyses like the one in our case study are important for defining which groups could benefit the most from interventions, education, and policy changes when attempting to mitigate public health challenges. You can see in this article that many additional considerations would be involved to adequately understand the data enough to recommend policy changes.

Suggested Homework

Students can predict air pollution monitor values using a different algorithm and provide an explanation for how that algorithm works and why it may be a good choice for modeling this data.

Helpful Links

1. A review of tidymodels
   
2. A course on tidymodels by Julia Silge
   
3. More examples, explanations, and info about tidymodels development from the developers
   
4. A guide for pre-processing with recipes
   
5. A guide for using GGally to create correlation plots
   
6. A guide for using parsnip to try different algorithms or engines
   
7. A list of recipe functions
   
8. A great blog post about cross validation
   
9. A discussion about evaluating model performance for a deeper explanation about how to evaluate model performance
   
10. RStudio cheatsheets
11. An explanation of supervised vs unsupervised machine learning and bias-variance tradeoff.
12. A thorough explanation of principal component analysis.
13. If you have access, this is a great discussion about the difference between independence, orthogonality, and lack of correlation.
14. Great video explanation of PCA.

Terms and concepts covered:
Tidyverse
parameters and hyperparameters

Packages used in this case study:

Package	Use